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Preface 


The  purpose  of  this  study  was  to  investigate  the  nonlinear  behavior  of  a 
simple  po'-ared  1  .'.sting  hypersonic  vehicle  flying  in  a  near  circular  orbit 
above  a  spherical  nonrotating  Earth  with  gradients  in  atmospheric  density 
and  pressure  and  an  inverse  square  law  for  gravity.  The  vehicle  is 
constrained  to  fly  in  a  vertical  plane  so  only  longitudinal  motion  is 
modeled.  Bifurcation  analysis,  utilizing  the  AUTO  software  package,  was 
used  to  conduct  this  study.  A  simple  five-state  model  with  three 
different  thrust  laws  was  used  to  describe  an  unaugmented  vehicle  whose 
geometric  and  aerodynamic  characteristics  follow  those  of  the  literature. 
A  parameter  represented  a  body  flap  deflection  (5^)  was  used  to  conduct 
one  set  of  bifurcation  sweeps  for  each  thrust  law.  Then  a  second  set  of 
bifurcation  sweeps  for  each  thrust  law  was  obtained  using  a  parameter 
representing  a  throttle  (6T)  which  scaled  the  value  of  the  thrust. 
Secondary  parameters  representing  simple  feedback  gains,  were  subsequently 
added. 

I  wish  to  extend  my  sincerest  thanks  to  my  thesis  advisor,  Capt  Jim 
Planeaux,  for  his  patient  and  caring  nature.  His  guidance  and  insight 
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while  reviewing  this  document.  Finally,  but  most  importantly,  I  would  like 
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a  :  angle  of  attack  (radians) 

8bf  :  body  flap  deflection  (degrees) 
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p  :  longitude  (radians) ,  also  eigenvalue 

6  r  pitch  angle  (radians) 
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:  angular  velocity  of  body  frame  to  inertial  frame  (radians/sec) 

:  angular  velocity  of  Earth  (radians/sec) 
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Cpo  :  lift  coefficient  at  zero  angle  of  attack 

Cpg  :  partial  derivative  of  lift  w.r.t.  angle  of  attack  (1/radian) 

C,  :  pitching  moment  coefficient 
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:  partial  derivative  of  pitching  moment  w.r.t. 

body  flap  deflection  (1/deg) 

:  partial  derivative  of  pitching  moment  w.r.t. 

angle  of  attack  (1/radian) 

:  partial  derivative  of  pitching  moment  w.r.t. 

pitch  rate  (sec/radian) 

:  drag  (lb) 

:  force  generated  by  constant  torque  (F=T/r)  in  example  problem  (lb) 
:  gravity  as  a  function  of  radius  from  Earth's  center  (ft  -  sec’^) 

:  altitude  (ft) 

:  (I,  -  I,) /I, 

:  radius  of  gyration  in  pitch  (ft) 

:  pitch  rate  feedback  gain 
:  characteristic  length  (vehicle  length)  (ft) 

:  lift  (lb) 

:  mass  (slugs) 

:  sum  of  moments  about  vehicles  center  of  mass  (ft  -  lb) 

:  Mach  number 
:  ambient  pressure  (psf) 

:  nozzle  exit  pressure  (psf) 

:  pitch  rate  of  vehicle  relative  to  Earth  (radian/sec) 

:  pitch  rate  relative  to  the  Earth  at  starting  altitude  (radian/sec) 
:  geocentric  radius  (ft),  also  length  of  rod  in  example  problem  (ft) 
:  area  of  lifting  surface  (ft*) 

:  thrust  (lb),  Period  (seconds),  torque  (ft-lb) 

:  velocity  (ft/sec) 

:  rocket  nozzle  exhaust  velocity  (ft/sec) 

:  weight  (lb) 
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Xr  :  partial  derivative  of  the  longitudinal  force  w.r.t.  radius 

Xu  :  partial  derivative  of  the  longitudinal  forces  w.r.t.  velocity 

in  the  &x  direction 

()e  :  equilibriuni  value;  also  value  for  example  problem  (i.e.  ©e) 

{)0  ;  values  at  initial  starting  equilibrium  solution 

(°)  :  derivative  with  respect  to  time 
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Abstract 

Bifurcation  analysis  was  used  to  investigate  the  nonlinear  behavior 
of  a  simple  powered  lifting  hypersonic  vehicle  in  circular  orbit  about  a 
spherical  nonrotating  Earth  with  gradients  in  atmospheric  density  and 
pressure  and  an  inverse  square  law  for  gravity.  Vehicle  motion  is 
constrained  to  a  vertical  plane  so  only  longitudinal  dynamics  were 
modeled.  Bifurcation  analysis  was  conducted  using  the  AUTO  software 
package.  A  simple  five-state  model  with  three  different  thrust  laws  was 
derived  to  describe  an  unaugraented  vehicle  whose  geometric  and  aerodynamic 
characteristics  follow  those  of  the  literature.  A  parameter  representing 
a  body  flap  deflection  {5bj)  was  used  to  conduct  one  set  of  bifurcation 
sweeps  for  each  thrust  law.  A  second  set  of  bifurcation  sweeps  for  each 
thrust  law  was  obtained  using  a  parameter  representing  a  throttle  (5T) 
which  scaled  the  thrust.  Secondary  parameters  representing  simple 
feedback  gains  were  subsequently  added.  Results  were  surprising  for  a 
simple  system  with  basically  linear  aerodynamics.  Periodic  branches 
arising  from  the  loss  of  pitch  stability  or  associated  with  a  "resonance 
altitude"  are  routinely  found  with  significant  amplitude,  and  periods  on 
the  order  of  an  elliptical  orbit's  period  for  a  given  geocentric  radius. 
Rotational  states  generally  had  sub-oscillations  of  greater  frequency. 
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BIFURCATION  ANALYSIS  OF  THE  LONGITUDINAL  DYNAMICS 
OF  A  SIMPLE  POWERED  LIFTING  HYPERSONIC  VEHICLE 


I.  Introduction 


Introductory  Discussion 

In  the  past  interest  in  hypersonic  vehicle  dynamics  has  concentrated 
on  assuring  stable  reentry  and  return  to  a  few  specific  points  on  the 
earth.  The  need  for  maneuverability  was  limited.  The  renewed  interest  in 
lifting  hypersonic  vehicle  dynamics  and  design;  brought  about  by  the 
Trans-Atmospheric  Vehicle  projects,  Boost  Glide  Vehicles,  the  National 
Aerospace  Plane  and  other  hypersonic  lifting  vehicles  designed  for 
improved  maneuverability  and  greater  versatility;  has  generated  a  need  to 
better  understand  and  anticipate  their  possible  nonlinear  dynamic  effects. 
Specifically,  the  ability  to  predict  nonlinear  dynamic  responses,  periodic 
equilibrium  states  and  other  dynamic  phenomena  for  a  representative 
hypersonic  vehicle  is  growing  in  importance.  Previous  work  by  Etkin  (4), 
Berry  (2),  Vihn  (15)  and  others  demonstrate  some  interesting  behavior  of 
the  longitudinal  dynamics  for  hypersonic  vehicles  due  to  the  variation  of 
atmospheric  density,  gravity  and  Mach  number  with  altitude.  Bifurcation 
and  continuation  analyses  have  been  used  successfully  to  examine  the 
nonlinear  behavior  of  fighter  aircraft  in  a  variety  of  configurations.  It 
was  felt  this  type  of  analysis  would  yield  insight  into  the  nonlinear 
behavior  of  a  hypersonic  vehicle  as  well. 
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The  purpose  of  this  thesis  was  to  explore  the  nonlinear  dynamic 
responses  of  a  hypersonic  vehicle  and  using  a  more  global  technique  to 
investigate  these  effects.  In  addition,  it  is  hoped  the  application  of 
bifurcation  analysis  techniques  to  the  highly  nonlinear  hypersonic  regime 
would  help  extend  the  basic  techniques  available  for  further  analysis  of 
hypersonic  vehicles. 

Summary  of  Previous  Studies 

Several  papers  have  been  presented  over  the  last  four  decades  that 
impact  directly  on  the  study  of  the  longitudinal  dynamics  of  a  hypersonic 
vehicle  flying  in  an  atmosphere  that  contains  gradients  in  density  and  the 
effects  of  curvature  of  the  flight  path.  Host  of  these  previous  works 
built  in  some  way  upon  the  work  presented  in  1950  by  Neumark  (9),  which 
was  then  extended  to  a  lifting  vehicle  in  orbital  flight  by  Etkin  (4)  in 
1961.  The  results  of  these  later  works  have  served  to  enhance  the 
material  originally  found  in  these  two  landmark  papers  for  lifting 
vehicles.  Having  said  this  one  should  note  that  some  correction  of 
Etkin ' s  observations  regarding  the  behavior  of  the  phugoid  and  pitching 
mode  characteristics  are  found  in  the  work  by  Vihn  and  Dobrzelecki  (15) 
and  verified  in  a  later  paper  by  Harkopoulos,  et  al  (7)  as  well  as  this 
author's  most  recent  work. 

In  his  paper  Neumark  details  the  motivation  for  bis  work  which  was 
based  on  several  of  the  very  first  studies  of  the  behavior  of  airplanes  in 
steep  angled  dives.  He  is  one  of  the  original  writers  on  the  subject  of 
the  effect  of  density  gradients  on  airplanes  having  published  his  first 
work  on  the  subject  in  1931;  the  earliest  being  in  1929.  His  paper 
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published  in  1950  was  the  first  published  work  in  which  the  longitudinal 
equations  of  motion  were  cast  in  the  form  of  a  quintic?  having  added  an 
equation  to  describe  the  change  in  altitude  with  time.  This  form  of  the 
longitudinal  equations  gives  rise  to  a  fifth,  real  root  (eigenvalue)  which 
gives  an  indication  of  the  vehicles  ability  to  hold  a  fixed  equilibrium 
altitude.  The  results  of  his  study  demonstrate  the  increasing  affect  the 
density  gradient  has  on  the  longitudinal  dynamics  as  the  speed  of  the 
vehicle  approaches  Mach  one.  Neumark  found  that  the  principal  effect  of 
the  density  gradient  is  on  the  phugoid  mode?  he  states  that  the  short 
period  (pitching)  mode  is  unaffected.  Increasing  density,  increases  the 
phugoid  frequency  thus  shortening  the  period.  The  effect  on  phugoid 
damping  was  not  clear  having  been  complicated  by  compressibility  effects 
at  speeds  above  H  *1.4  .  He  concluded  the  height  mode  would  have  a  very 
long  time  constant  and  may  be  either  a  subsidence  or  divergence  and  has 
importance  only  for  hypersonic  flight  or  flight  at  constant  altitude  for 
long  periods  of  time  (9:325). 

Etkin's  classic  of  1961  extends  Neumark's  work  to  the  truly 
hypersonic  case  and  includes,  necessarily,  the  mathematical  modifications 
to  account  for  the  curvature  of  the  undisturbed  flight  path  and  the 
variation  of  gravity  with  altitude.  In  his  analysis  the  longitudinal 
equations  of  motion  for  flight  in  a  vertical  plane  about  a  nonrotating 
Earth  whose  atmosphere  is  at  rest  are  linearized  and  the  behavior  of  the 
vehicle  subject  to  small  disturbances  about  an  equilibrium  is  examined. 
In  addition,  he  presents  results  from  numerical  solutions  to  the  nonlinear 
equations  and  does  a  comparison  with  the  linear  approximation.  Of  note  in 
the  equations  of  motion  is  the  addition  of  the  torque  about  the  vehicle's 
center  of  mass  due  to  the  small  variation  of  gravity  acting  on  a  body  at 


3 


very  high  altitudes  (above  500,000  ft)  where  the  pitch  damping  is 
negligible.  In  this  realm  the  gravity  torque  generates  the  dominant 
moment  acting  on  the  vehicle  and  for  a  standard  vehicle  configuration 
whose  longitudinal  axis  is  nearly  aligned  with  the  flight  path  this  effect 
is  destabilizing. 

Etkin  examines  four  basic  cases  using  the  same  (steady  reference) 

lift  coefficient  (Ct  =  0.05  [rad'*]). 

Case  A  Constant  thrust  rocket,  full  set  of  equations 
Case  B  Air-breathing  engine  (T  «  p),  full  set  of  equations 

Case  C  Approximate  equations  (i.e.  no  density  gradient) 

Case  D  Constant  thrust  rocket  with  qj  =  0 

where  qg  is  the  equilibrium  value  of  the  pitch  rate  relative  to  the  Earth. 

Etkin  determined  that  the  effects  of  varying  density  and  gravity 

with  altitude  and  the  effects  due  to  the  Earth's  curvature  and  the  thrust 

law  have  significant  impact  on  the  phugoid  mode  and  the  stability  of  the 

height  mode,  but  insignificant  effect  on  the  pitching  (short  period)  mode 

except  at  very  high  altitude  where  the  pitch  damping  is  overcome  by  the 

gravity  torque.  In  addition,  Etkin  found  that  above  400,000  ft  the  period 

of  the  pitching  and  phugoid  modes  approached  each  other  and  he  asserts 

that  they  become  equal,  after  which  the  phugoid  tends  toward  the  orbital 

period  and  the  pitching  mode  tend  toward  infinity.  In  this  altitude  range 

he  demonstrated  a  dynamic  coupling  between  the  two  modes  and  when  nearly 

equal  all  relation  to  two  classical  modes  breaks  down  with  substantial 

pitching  motion  in  the  phugoid  mode.  Finally  when  the  two  modes  are 

exactly  equal  he  determines  the  system  to  be  unstable.  For  the  height 

mode  Etkin  found  that  it  represents  "a  spiral,  proceeding  away  from  the 

reference  orbit."  He  also  noted  an  interesting  variation  in  the  way  the 

speed  changes  with  altitude  above  and  below  350,000  ft.  Above  this 
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altitude  as  the  altitude  is  decreased  the  speed  increases  whereas  below 
this  altitude  the  opposite  is  the  case  (4:787-738). 

In  the  work  by  Vihn  and  Dobrzelecki  (15)  an  "analytic  study  of  the 
longitudinal  dynamics  of  a  thrusting,  lifting  orbital  vehicle  in  a  nearly 
circular  orbit"  is  conducted.  The  basic  set  of  five  equations  used  to 
describe  the  longitudinal  motion  of  a  vehicle  in  orbit  about  a  spherical 
Earth  were  used.  A  strictly  linear  analysis  as  well  as  analysis  including 
second  order  terms  in  the  Taylor  series  expansion  of  the  atmospheric 
density  were  used  to  develop  explicit  relationships  to  describe  the 
orbital  motion.  Also  developed  were  analytic  expressions  for  the  period 
and  damping  of  the  "angle  of  attack"  (pitching)  mode.  As  with  Etkin  they 
observed  an  altitude  where  the  velocity-altitude  relationship  inverts. 
They  went  on  to  develop  an  expression  based  on  vehicle  characteristics 
that  defines  the  altitude  where  this  "inversion"  takes  place.  Finally 
they  found  the  trend  at  high  altitude  of  the  linearized  phugoid  or  long- 
period  mode  and  angle  of  attack  (pitching)  mode  tend  to  become  nearly 
equal  in  frequency,  period  and  damping,  then  diverge.  (Figure  1)  Similar 
behavior  for  the  very  same  equations  was  found  earlier  by  Etkin  (4:785- 
788)  where  he  concluded  the  two  modes  "crossover"  and  the  phugoid  period 
tended  to  the  orbital  period  while  the  period  of  the  pitching  (short 
period)  mode  tended  toward  infinity.  At  the  point  of  "crossing"  Etkin 
concluded  the  dynamic  system  would  be  unstable  (4:787). 

Stengel  also  found  the  same  basic  trends  in  the  three  longitudinal 
modes  however  his  work  looked  more  closely  at  the  stability  questions  and 
dealt  at  some  length  with  various  techniques  to  provide  altitude  stability 
for  a  vehicle  in  supersonic  cruising  flight  (13).  In  his  work  he  uses  the 
linearized  equations  for  longitudinal  motion,  characterized  by  the 
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Figure  1.  Frequencies  near  Resonance 
Altitude  (15:19) 

perturbation  variables  Au  (forward  velocity),  Aa  (angle  of  attack),  A© 
(pitch  attitude) ,  and  Ah  (altitude) ,  to  study  the  interrelations  of 
motions  these  variable  characterize.  He  then  used  this  information  to 
test  how  various  combinations  of  feedback  and  control  would  affect 
altitude  stability.  In  addition  to  developing  analytical  transfer 
functions  he  conducted  numerical  studies  and  summarizes  the  effectiveness 
of  the  various  techniques  proposed.  Some  of  his  results  were  tabulated  in 
his  work  and  are  reproduced  on  the  following  page  in  Table  1. 

In  a  later  work  by  Berry  (2),  he  examined  the  effect  on  the  "long- 
period  dynamics"  of  a  vehicle  of  similar  characteristics  to  previous 
authors  but  included  an  "advanced  air-breather"  with  a  more  complicated 
thrust  law  which  is  more  representative  of  a  true  hypersonic  vehicle  like 
the  National  Aerospace  Plane.  He  concluded  that  the  height  mode  stability 
and  long-period  damping  were  strongly  affected  by  the  slope  of  thrust  with 
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Mach  number  (6:255-257).  While  all  the  trends  indicated  in  his  work  are 
valid  they  are  not  in  the  strictest  sense  complete.  It  is  the  variation 
of  the  difference  between  thrust  and  drag  with  height  and  Mach  number  that 
determines  stability.  These  more  complete  relations  are  identified  in 
works  by  Markopoulos,  et  al  (7:285)  and  Stengel  (13:468,472).  Berry 
further  examined  the  effectiveness  of  various  simple  feedback  schemes 
involving  only  the  pitch  control  surface  to  stabilize  the  long-period 
and/or  height  modes.  Some  results  for  feeding  back  combinations  of  pitch 
attitude,  forward  velocity  and  altitude  to  the  pitch  control  surface  for 
the  rocket  and  "advanced  air-breather”  were  presented  (2:256-257). 

Most  recent  is  the  work  by  Markopoulos,  Mease  and  Vihn  (7)  where  the 
linearized  equations  of  motion  are  used  to  examire  the  thrust  law  effects 
on  the  longitudinal  dynamics  of  an  aerospace  vehicle  flying  at  hypersonic 
speeds.  Their  work  demonstrates  the  dependence  of  the  height  mode 
stability  and  phugoid  damping  on  the  way  the  longitudinal  force  varies 
with  both  altitude  and  speed.  In  addition  they  confirm  the  results  of  the 
previous  study  by  Vihn  and  Dobrzelecki  (15)  for  the  high  altitude  trend  of 
the  phugoid  or  long-period  mode  and  angle  of  attack  (pitching)  mode. 

Of  special  interest  in  the  work  by  Markopoulos,  et  al,  is  the 
characterization  of  expected  height  mode  stability  and  phugoid  damping 
over  the  plane  of  all  thrust  possibilities.  The  plane  is  defined  by  two 
parameters,  specifically  the  variation  in  the  longitudinal  force  with 
height  (Xr)  and  velocity  (Xu) .  The  correlation  of  the  points  relating  to 
earlier  work  by  Etkin  (4)  are  in  excellent  agreement.  They  went  on  to 
conclude  that  it  is  actually  "the  partial  derivatives  of  the  difference 
between  thrust  and  drag  with  respect  to  speed  and  altitude  that  plays  the 
key  role  in  determining  the  stability  of  the  translational  dynamics  (7:287)." 
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TABLE  1 


SOME  FEEDBACK  EFFECTS  FOR  STENGEL'S  STANDARD  CASE  WITH  AUGMENTED  SHORT 
PERIOD  (FEEDBACK  IS  NEGATIVE  UNLESS  DENOTED  BY  (+)) 


Feedback  Variable 


and  Control 

Height  Mode 

Phugoid  Mode 

Attitude  to: 

thrust 

SS 

SI 

lift  (-) 

I 

S 

lift  (+) 

SS 

I 

moment  (-) 

S 

I 

moment  (+) 

SI 

S 

Pitch  angle  to: 

thrust 

N 

I 

lift  (-) 

I 

SS 

lift  (+) 

N 

SI 

moment  (-) 

I 

ssa 

moment  (+) 

N 

SI 

Forward  Velocity  to: 

thrust 

SS 

ssa 

lift  (-) 

SI 

SS 

lift  (+) 

N 

SI 

moment 

SS 

ssa 

Angle  of  attack  to: 

thrust 

N 

SS 

lift  (-) 

SI 

SIa 

lift  (+) 

N 

SS 

moment  (-) 

? 

SSb 

moment  (+) 

? 

S  =  Stability  SS  =  Strong  Stability 

I  =  Instability  SI  =  Strong  instability 
N  =  Neutral  stability 
3  With  Limited  Travel. 

b  Conditional  Stability  Stengel  (13:470) 

mmmmammmmmmmmmmmmmmmmmmmmmmmmmmmmmaMummmmmmmmmmamamBmmmmHmmmmmmmmmmmmam 


This  same  relationship  was  discussed  by  Stengel  in  his  article  on 
"Altitude  Stability  for  Supersonic  Cruising  Flight"  (13:468,472). 
Markopoulos,  et  al,  further  concluded  after  numerical  simulation  of  their 
full  and  reduced  order  mathematical  models  that  over  the  plane  of  all 
engine  possibilities 

"if  thrust  increases  faster  than  drag  vith  respect  to  speed  at 
least  one  of  the  translational  modes  (height  or  phugoid)  vill 
be  unstable.  Increasing  the  partial  derivative  of  the  differ- 


8 


ej ice  between  thrust  and  drag  with  respect  to  altitude  has  a 
destabilizing  effect  on  the  height  mode  and  a  stabilizing 
effect  on  the  phugoid  (7:287)." 

Finally,  they  concluded  that  to  first  order,  the  period  of  the  phugoid 
mode  as  well  as  all  characteristics  of  the  pitching  mode  are  independent 
of  the  thrust  law  (7:287). 

Outline  of  Analysis 

In  this  thesis  the  dynamic  behavior  of  a  powered  lifting  hypersonic 
vehicle  in  nearly  circular  orbital  flight  about  the  center  of  mass  of  a 
spherical  nonrotating  planet  (specifically  the  Earth)  whose  atmosphere 
contains  gradients  in  density  and  pressure  and  whose  gravitational  field 
follows  the  inverse  square  law  is  examined. 

Throughout  this  thesis  emphasis  is  given  to  the  leading  order 
aerodynamic  behavior  and  simplifying  assumptions  to  this  end  are  brought 
to  the  readers  attention  as  required.  In  order  to  focus  the  scope  of  this 
work  it  is  assumed  the  flow  is  inviscid  therefore  the  effects  of  high 
temperature  gas  flows  are  neglected.  This  assumption  is  consistent  with 
general  longitudinal  analysis  found  routinely  in  the  literature  and  allows 
use  of  simple  Newtonian  impact  theory  as  the  basis  for  the  aerodynamics. 

To  begin  this  study  the  reader  should  have  a  good  mental  image  of 
the  problem  being  analyzed,  and  a  well  developed  understanding  of  the 
equations  used  to  describe  the  translation  and  rotation  of  a  body  flying 
a  great  circle  about  a  spherical  planet.  To  facilitate  this  the  basic 
equations  for  a  vehicle  flying  in  an  atmosphere  at  rest  relative  to  a 
nonrotating  spherical  planet  (tt?=0)  are  derived.  The  first  step  in 
analyzing  this  problem  is  to  identify  an  inertial  reference  frame.  Then, 
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three  advantageous  frames  of  reference  relative  to  the  inertial  frame  are 
introduced  from  which  a  set  of  equations  describing  the  forces  and  moments 
acting  on  the  body  of  interest  are  developed.  It  is  common  in  trajectory 
analysis  to  use  a  wind  axis  system  as  shown  in  Figure  2  where  the  positive 
x-axis  is  parallel  to,  but  opposite,  the  relative  wind.  In  this  way  the 
aerodynamic  forces  are  cleanly  defined  and  the  velocity  vector  has  a 
single  non-zero  component.  For  the  analysis  of  angular  momentum  the  body- 
fixed  axes  are  used,  also  shown  in  Figure  2  as  bx  and  bz ,  thus  the  moments 
of  inertia  are  time  invariant  and  for  a  fixed  mass  and  mass  distribution, 
as  is  the  case  here,  the  moments  of.  inertia  are  constant.  The  vehicle 
axes  indicated  in  Figure  2  by  (v)_  ,  is  used  as  a  convenient  intermediate 
frame  between  the  body  or  wind  axes  and  the  inertial  axes. 

As  with  linear  analysis,  the  bifurcation  analysis  begins  at  a  known 
equilibrium  point,  but  rather  than  linearizing  the  equations  and  looking 
at  small  disturbances  about  this  point,  a  continuation  method  is  used  to 
solve  for  the  flow  of  equilibrium  solutions  (specifically  the  pseudo-arc 
length  technique  resident  to  AUTO;  the  software  package  used  for 
continuation  and  bifurcation  analysis  in  this  study)  (3:12-16;  12:116). 
From  the  path  of  equilibrium  solutions  (or  stationary  points)  bifurcating 
solutions,  or  in  other  words,  additional  solution  paths  are  located  and 
explored.  Of  the  various  types  of  possible  solution  branches  special 
interest  is  given  to  branches  obtained  subsequent  to  a  nonhyperbolic  or 
degenerate  point  (12:18-20).  These  often  give  rise  to  branches  of 
periodic  solutions  where  motion,  such  as  periodic  oscillations  develop. 
On  these  branches  the  dynamic  behavior  comes  to  life.  Many  of  these 
concepts  are  clarified  in  section  II  where  the  nature  of  bifurcation 
analysis  is  discussed. 
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Equations  of  Motion 


The  author  is  indebted  to  Planeaux  (10) ,  HcRuer,  et  al  (8:204,220) 
and  Etkin  (5:104,148)  for  background  and  guidance  for  the  following 
development.  The  reader  is  encouraged  to  review  the  two  texts  for  details 
on  the  development  of  the  following  system  of  equations.  To  begin,  the 
basic  assumptions  on  which  the  equations  of  motion  are  based  must  be 
stated: 

1.  The  Earth  is  an  inertial  reference  frame. 

2.  The  vehicle  is  a  rigid  body. 

3.  The  vehicle  mass  and  mass  distribution  are  constant. 

4.  The  vehicle  is  symmetric  about  the  x-z  plane. 

5.  The  body  fixed  axes  are  aligned  with  the  principal  axis  of  the 
vehicle. 

In  addition  to  the  simplifications  resulting  from  the  assumptions 
above,  the  terms  associated  with  motion  in  the  horizontal  plane  are 
neglected  leaving  the  system  as  shown  above  in  Figure  2. 

Looking  first  at  the  angular  and  kinematic  relations,  by  inspection  the 
following  angular  rate  of  the  three  reference  frames  relative  to  the 
inertial  frame  are  given  as: 

=  gj  =  (-£  +  $)  j  »  (6  +  -£)  £y 


=  (-£)  J  =  "|i  tfy 


(2) 
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(3) 


s±wi  =  (-£  +  Y)  J  =  (Y  -  A)  ^y 


where:  (i,  j,£) 

(£xlSyl£s) 

(tix,fiyl&2) 

Y 

6 


i*vi 

q 


unit  vectors  of  inertial  frame 

unit  vectors  of  body  frame 

unit  vectors  of  wind  frame 

time  rate  change  of  longitude  (rad/sec) 

time  rate  change  of  flight  path  angle 

(rad/sec) 

time  rate  change  of  pitch  angle  (rad/sec) 
angular  velocity  of  body  to  inertial  frame 
(rad/sec) 

angular  velocity  of  vehicle  to  inertial 
frame  (rad/sec) 

angular  velocity  of  wind  to  inertial  frame 
(rad/sec) 

pitch  rate  of  the  vehicle  relative  to  the 
Earth  (rad/sec) 


The  radius  (r)  is  the  distance  from  the  center  of  mass  of  the  Earth  to  the 
vehicles  center  of  mass  and  written  as  a  vector  in  the  vehicle  frame  is 
given  as: 


r  = 


(4) 


where:  =  unit  vectors  of  the  vehicle  frame 

r  =  geocentric  radius  (ft) 

From  eqn(4)  the  velocity  can  be  written  by  differentiating  in  the  vehicle 

frame: 


V  =  £  =  -£VZ  +  iiivixr 
=  -i$2  + 


where:  y 

V 


flight  path  angle  (rad) 
velocity  (ft/sec) 
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However  the  velocity  can  also  be  shown  to  be: 

V=  V(cos  (y)  Vx  -  sin(y)  Qg)  (6) 

By  combining  eqn(5)  and  eqn(6)  the  following  scalar  kinematic  relations 
are  obtained: 

i  =  Vsin(y)  (7) 

|i  =-pCos(y)  (8) 

where:  y  =  flight  path  angle  (rad) 

r  =  radius  (ft) 

V  =  velocity  (ft/sec) 

From  eqn(l)  the  equation  for  the  time  rate  change  of  pitch  angle  can  be 
obtained  by  solving  for  6  as  follows: 

6  =  g  +  |i  (9) 

Following  directly  from  the  assumptions  3  through  5  the  time  rate  change 
of  the  pitch  velocity  of  the  vehicle  is  given  as  (4:145): 

$  =  -y-  (10) 

where:  Iy  =  moment  of  inertia  about  the  y-axis  of  the  body  (slug-ftJ) 

=  sum  of  moments  about  the  vehicle^  center  of  mass  (ft-lb) 

Equations  (7,9  and  10)  make  up  three  of  the  five  equations  of 

motion.  The  remaining  two  equations  fall  out  of  the  force  balance  which 

is  dealt  with  next. 
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Summing  forces  acting  at  the  vehicles  center  of  mass  yields: 


F  -  m&  =  mg(z)  Vz  +  tSx  -  L$z  -  D$x 

=  [-mg(r)  sin(y)  +  Tcos(a)  -  D]#x 

+  Dn£rcos(y)  -  rsin(o)  -  L] 


(11) 


where:  a  =  acceleration  <ft-sec"2) 

D  =  drag  (lb) 

g(r)  =  gravity  as  a  function  of  geocentric  radius  (ft-sec-2) 

L  =  lift  (lb) 

m  =  mass  (slug) 

r  =  geocentric  radius  (ft) 

T  =  thrust  (lb) 

Now  velocity  in  the  wind  frame  is  written  as  V=Vtix.  The  acceleration 
in  the  wind  frame  is  given  by: 


V  =  a  =  V#x  -  (y-(x)  Vtiz 


(12) 


Since  F=ma  eqn(12)  and  eqn(ll)  can  be  set  equal,  after  multiplying 
eqn(12)  by  the  mass  m,  and  upon  separating  into  scalar  components  yields 
the  following  two  equations. 


mV  =  Tcos(a)  -  D  -  mgsin(y) 


(13) 


-mtf  -  |i)  V  -  -Tsin(a)  -  L  +  mg cos(y)  (14) 


Utilizing  the  relation  0  =  y  +  a  and  eqn(9)  the  following  expression  is 
obtained: 


g-d  = 


(15) 
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Substituting  eqn(15)  into  the  left  hand  side  of  eqn(14)  and  solving  ford 
yields  the  final  equation  for  the  set  of  five  equations  of  motion. 

&  =  g  -  -4sin(«)  -  -4  +  J^ilicostY)  (16) 

mV  mV  r 

Equations  (7,9,10,13  and  16)  comprise  the  set  of  five  dynamic  and 
kinematic  equations  for  this  analysis.  Together  with  the  following 
expressions  for  the  aerodynamic  forces  and  moments,  and  the  thrust 
equations  they  comprise  the  complete  set  of  equations  required  to  conduct 
this  study.  The  five  equations  of  motion  are  reprinted  below  for 
convenience: 


mV  =  Tcos(a)  -  D  -  mgs in(y) 

(13) 

&  =  G  -  -^sin(a)  -  —  +  “4-cos (y) 
mV  mV  i 

(16) 

6  =  g  +  A 

(9) 

•Ly 

(10) 

£  =  Vsin(y) 

(7) 

Aerodynamic  Forces  and  Moment  Coefficients.  As  with  linear  analysis 
the  standard  forms  cf  the  forces  and  the  moment  due  to  aerodynamics  will 
be  used,  and  are  listed  below.  Notice  however,  the  term  on  the  right  hand 
side  of  eqn(19).  This  term  is  the  moment  about  the  center  of  mass  of  a 
satellite  in  a  gravitational  field  and  as  found  by  Etkin  is  a  significant 
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factor  at  high  altitudes  (4:783;  11:21) 


D  = 

(17) 

L  = 

--  f  P  V“SCt 

(18) 

*-■§ 

P  v2sicb 

-  A  £L  ( j  -  j  ) 
2  r  x  3 

sin  (20) 

(19) 

ii 

Cp0  +  9D6 

(20) 

CL  * 

C LO  +  C La  a 

(21) 

^JD  =  CjO 0 

+ 

*  (g  -  ji)  + 

CMribf 

(22) 

where: 


p  =  atmospheric  density  (slug/ft3) 

CD  =  nondimensional  drag  coefficient 

CD0  =  basic  drag  coefficient 

CDa  =  partial  derivative  of  Cp  w.r.t.  alpha  (1/rad) 

=  nondimensional  lift  coefficient 
C^q  =  basic  lift  coefficient 

CLa  =  partial  derivative  of  Ct  w.r.t.  alpha  (1/rad) 

Cj  =  nondimensional  aerodynamic  moment  coefficient 

Cj,  =  partial  derivative  of  C|  w.r.t.  pitch  rate  (sec/rad) 

Cjg  =  partial  derivative  of  Ca  w.r.t.  alpha  (1/rad) 

C,tbf  =  partial  derivative  of  C.  w.r.t.  body  flap  deflection 

(1/deg) 

g(r)  =  gravitational  acceleration  (ft/sec2) 

1  =  characteristic  length  (vehicle  length)  (ft) 

S  =  area  of  lifting  surface  (ft2) 


The  final  term  in  eqn(22)  represents  the  contribution  of  the  body  flap 
deflection  to  the  moment  coefficient  and  is  used  as  a  standard  pitch 


control  surface.  The  density  is  calculated  using  one  of  two  analytic 
expressions  depending  on  the  altitude.  The  specifics  of  how  the  density 
is  calculated  as  well  as  a  brief  discussion  of  the  development  of  the 


analytic  expressions  is  found  in  Appendix  2. 
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Thrust  Laws.  Three  basic  ideal  thrust  laws  are  used  in  this  study. 
The  equations  representing  these  thrust  laws  are  detailed  below. 

constant  thrust  rocket 

T  =  -|p0  vfcCi*  (23) 

where:  pc  =  atmospheric  density  at  starting  altitude  (slug/ft3) 

Vq  =  velocity  at  starting  altitude  (ft/sec) 

Cj,q  =  drag  coefficient  at  starting  altitude 

Notice  for  the  constant  thrust  case  the  thrust  is  fixed  at  the  values  of 

the  drag  for  the  starting  altitude  (i.e.  the  altitude  where  the 

bifurcation  sweeps  starts).  This  is  required  since  to  start  the 

continuation  method  an  equilibrium  solution  must  be  provided  as  a  first 

step.  In  all  cases  the  equilibrium  solution  has  a  =  9  =  0  radians.  This 

requires  the  thrust  to  equal  the  drag  for  equilibrium. 

variable  thrust  rocket  (6:356) 

T  =  itiV^  +  <P0  -  Pa)  Ad  (24) 

As  stated  above  the  starting  equilibrium  point  with  a  =  0  =  0  radians 
requires  the  thrust  to  equal  drag  when  the  continuation  method  is  begun. 
This  requires  that  the  mass  flow  be  determined  by  setting  the  thrust  equal 
to  the  drag  at  the  starting  altitude,  therefore  the  mass  flow  rate  of  the 
exhaust  is  fixed  at  the  following  equation.  Note  it  is  assumed  the  mass 
flow  is  sufficiently  small  relative  to  the  mass  of  the  vehicle  as  to  be 
negligible.  This  assumption  is  fairly  good  at  high  altitude  but  is  very 
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poor  at  altitudes  below  about  200,000  ft. 


ih  = 


-f  Po -  (P,  -  P.)  Ad 


(25) 


where: 


exhaust  nozzle  area  (ft^) 

mass  flow  rate  of  exhaust  (slug/sec) 

pressure  at  exhaust  nozzle  exit  plane  (psf) 

ambient  air  pressure  (psf) 

velocity  of  exhausted  mass  (ft/sec) 


ideal  turbojet  (9»332) 


T  = 


(26) 


As  with  the  constant  thrust  rocket,  thrust  for  the  ideal  turbojet  must 
equal  drag  at  the  starting  altitude  since  the  continuation  method  requires 
an  initial  "starting"  equilibrium  solution  and  a  =  8  =  0  was  taken  as 
the  values  of  these  states  at  the  equilibrium  point.  Therefore  the  value 
of  Tq  is  set  by  the  following  relation: 

T0  =  -|Po  VoSCoo  (27) 

Note  the  exponent  "x"  can  be  used  to  change  the  way  the  thrust  varies  with 
altitude.  For  the  standard  turbojet,  x  equals  one  (x=l). 

It  should  be  noted  that  in  all  cases  the  thrust  can  be  varied  with 
in  the  subroutine  CONST  through  a  Thrust  Scaling  Parameter  (5T)  that 
multiplies  the  calculated  value  of  thrust  using  the  above  relations. 
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This  parameter  thereby  acts  as  a  throttle,  increasing  or  decreasing  the 
thrust  as  required  for  the  ST  bifurcation  sweeps. 

Vehicle  Characteristics.  As  with  the  majority  of  previous  studies 
conducted  on  this  subject,  the  vehicle  geometry  and  aerodynamic 
characteristics  used  here  will  be  basically  the  same  as  those  used  by 
Etkin  (4:783-784).  This  allows  nearly  direct  comparison  of  results  which 
is  helpful  to  check  consistency  of  linearization  at  starting  points  and 
more  importantly  will  be  used  to  highlight  the  advantages  and  simplicity 
of  using  bifurcation  analysis  for  even  the  simplest  problems  in 
atmospheric  flight  mechanics/dynamics.  Etkin  obtained  his  data  from 
"simple  Newtonian  impact  theory  for  a  slender  body  (cone  or  wedge  of 
about  3°  semiangle)  at  moderate  angle  (4:784)."  In  this  study  a  small 
change  has  been  made  to  allow  for  lift  at  zero  angle  of  attack,  which  is 
more  representative  of  a  hypersonic  vehicle,  however  as  seen  in  Figure  3 
Etkin's  basic  lift  to  drag  ratio  was  followed  fairly  well. 

With  these  clarifications  stated  the  geometry  and  aerodynamic 
characteristics  are  as  follows: 


Geometry 

1  =  50  ft  k  =  (I  /m)1/2  =  25  ft  fca  (Ix -!,)/!,« -0.94 

W  =  700,000  lb  W7S  =  30  psf  (at  sea  level) 


Aerodynamics  (dimensions  are  [rad'*]  unless  otherwise  indicated) 


0.05 

0.50 

0.0133 

0.400 


'lO 

'»a 

;»q 

'i5bf 


0.00 

-0.0548 

-0.028 

-0.0822  [deg'1] 
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C  I  /  Cd 


Figure  3.  Comparison  of  Cp/Cp  with  Etkin's  work 
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II.  Introduction  to  Bifurcation  Analysis  and  Stability 


Equilibrium  Solutions  and  Equilibrium  Points 

At  this  point  it  is  best  to  discuss  some  of  the  central  points  of 
bifurcation  analysis  to  allow  the  uninitiated  to  understand  what  is  being 
done,  and  to  highlight  what  is  being  presented  in  the  bifurcation  diagrams 
and  phase  plane  diagrams  which  follow. 

Continuation  is  the  core  on  which  the  present  application  of 
bifurcation  analysis  is  based.  There  are  several  types  of  continuation 
techniques,  but  discussing  them  is  not  appropriate  here.  What  is  required 
is  a  general  description  of  what  is  obtained.  The  analysis  here  starts  at 
an  equilibrium  point  from  which  the  continuation  technique  is  begun.  The 
equilibrium  point  is  obtained  by  setting  the  time  derivatives  of  the  state 
variables  to  zero  thus  creating  a  set  of  homogeneous  equations  and  solving 
simultaneously  for  the  values  of  the  state  variables  that  satisfy  the 
homogeneous  equations.  Once  this  starting  point  is  obtained  the 
continuation  method  can  begin.  The  continuation  method  will  search  in  the 
vicinity  of  this  point  until  it  finds  another  point  which  satisfies  the 
set  of  homogeneous  equations.  This  process  continues  over  the  given  range 
of  the  specified  parameter  and  within  the  bounds  established  for  the 
states,  until  a  complete  parameter-dependent  family  of  equilibrium 
solutions  to  the  set  of  homogeneous  equations  has  been  found.  An 
important  feature  of  the  bifurcation  software  AUTO  is  the  ability  to  find 
the  equilibrium  solution  path  despite  running  into  singular  points  (limit 
points  and  bifurcation  points) .  AUTO  uses  a  pseudo-arclength  technique 
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which  allows  it  to  continue  working  despite  encountering  a  singular 
Jacobian  of  the  linearized  set  of  equations.  Other  less  robust  techniques 
fail  at  these  points  where  the  slope  of  the  solution  with  respect  to  the 
parameter  is  not  unique,  as  in  the  case  of  a  bifurcation,  or  undefined,  as 
with  a  limit  point.  The  reader  is  referred  to  the  user's  manual  for  AUTO 
and  the  text  by  Seydel  for  further  information  on  this  technique  (3:12-16; 
12:116+) . 

Simple  Nonlinear  Behavior: 

While  knowing  an  equilibrium  solution  branch  is  interesting,  the 
true  value  in  this  analysis  for  those  interested  in  nonlinear  effects  is 
the  accurate  location  of  limit  and  bifurcation  points.  The  reasons  for 
interest  in  these  singular  points  are  many.  In  the  case  of  a  bifurcation 
this  point  represents  the  intersection  of  other  solution  branch (es). 
Depending  on  the  type  of  bifurcating  point  there  is  the  potential  for 
complex  motion  arising  from  the  nonlinear  nature  of  the  problem.  In  the 
analysis  of  longitudinal  motion  of  a  powered  lifting  hypersonic  vehicle 
the  two  most  prevalent  singular  points  encountered  are  limit  points,  which 
may  give  rise  to  hysteresis  type  behavior  or  an  exchange  of  stability,  and 
Hopf  bifurcations,  which  generally  occur  in  this  study  when  the  phugoid 
mode  eigenvalues  cross  the  imaginary  axis  transversely  and  either  lose  or 
gain  stability.  Generally  for  the  analysis  of  the  longitudinal  dynamics 
of  a  powered  lifting  hypersonic  vehicle  the  Hopf  point  signals  the  loss  of 
stability  in  the  phugoid  mode.  A  Hopf  point  is  of  special  interest  in  the 
study  of  nonlinear  dynamics  as  the  behavior  subsequent  to  a  Hopf 
bifurcation  is  generally  characterized  by  increasingly  complex  periodic 
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motion  known  as  limit  cycles.  For  systems  with  three  or  more  degrees  of 
freedom  the  Hopf  bifurcation  may  also  be  the  first  step  in  the  direction 
of  chaos  (12:25) . 

While  the  mathematics  of  bifurcation  analysis  is  not  within  the 
scope  of  this  thesis  an  understanding  of  physical  processes  is.  To 
further  clarify  some  of  the  points  made,  and  to  provide  a  basis  for 
understanding  what  a  basic  bifurcation  study  is  all  about,  the  following 
simple  example  is  presented  (10) . 

Figure  4  shows  a  mass  (m)  connected  by  a  mass-less  rigid  rod  to  a 
mass-less  rigid  sleeve  which  is  around  a  spinning  shaft.  The  friction 
coefficient  between  the  spinning  rod  and  the  sleeve  is  constant  therefore 
a  constant  torque  is  generated  which  is  transmitted  to  the  mass  as  a  force 
(F)  via  the  mass-less  rigid  rod.  In-set  in  Figure  4  is  the  free-body 
diagram  for  the  mass  and  the  reference  axes.  Note  the  pendulum  is  held  at 
a  constant  angular  position  by  the  torque  applied  to  the  sleeve. 

The  scalar  equation  of  motion  for  the  mass  m  is: 

mib  =  F  -  mg sin(0)  f28) 


Where:  F 

g 

m 

r 

T 


force  generated  by  the  constant  torque  (F=T/r)  [lb] 

gravitational  acceleration  (ft/sec2) 

mass  (slugs) 

length  of  rod  (ft) 

torque  applied  to  sleeve  (ft-lb) 


At  a  point  of  equilibrium  the  left  hand  side  of  eqn  (28)  is  equal  to  zero, 
that  is  there  is  no  change  of  the  state  variable  when  the  forces  acting  on 
the  mass  are  in  equilibrium. 
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This  brings  out  an  important  point  to  remember  when  looking  at  bifurcation 
diagrams:  the  stationary  solution  path  is  made  up  of  equilibrium 
solutions  (equilibrium  points)  and  no  change  of  the  state  variables  is 
invol  ved. 

For  the  equilibrium  solutions  eqn  (28)  is  a  homogeneous  equation. 
In  the  case  of  the  powered  lifting  hypersonic  vehicle  there  is  a  set  of 
homogeneous  equations  that  are  solved  to  locate  the  equilibrium  solution. 
To  actually  conduct  the  bifurcation  analysis  a  parameter  must  be 
established  for  which  the  continuation  process  finds  a  solution  path.  For 
the  example  problem,  the  parameter  can  be  identified  as  A  =  F/(mg)  and 
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the  homogeneous  equation  can  then  be  put  in  the  form  of  eqn(29)  below: 


0e  =  sin'1(— )  (29) 

*  mg 

By  varying  A  from  -1  to  1,  a  plot  of  the  equilibrium  solutions  for  6e  is 
easily  obtained  even  with  a  hand  calculator.  The  bifurcation  analysis 
however,  yields  information  on  stability  by  linearizing  the  differential 
equation  (or  set  of  equations)  and  solving  for  the  eigenvalues  of  the 
subsequent  Jacobian  thus  providing  the  following  bifurcation  diagram. 

In  looking  at  the  bifurcation  diagram  it  is  seen  there  are  at  least 
two  equilibrium  solutions  for  each  A  in  the  open  set  (-1,1).  The  two 
points  corresponding  to  A=-l  or  1  are  called  limit  points.  It  is  clear  to 
see  that  limit  points  have  only  one  value  of  the  equilibrium  point  for  a 
given  parameter  and  are  points  where  the  equilibrium  solution  path  turns 
back  with  respect  to  the  parameter;  thus  the  slope  is  undefined.  Note 
also  that  to  one  side  of  a  limit  point  no  equilibrium  solutions  exist  yet 
on  the  other  side  two  equilibrium  solutions  exist  for  each  A.  In  the 
example  problem  the  limit  points  also  correspond  to  points  where  stability 
is  either  lost  or  gained,  depending  on  the  direction  of  the  parameter  A, 
however  this  is  not  always  the  case  for  limit  points  in  general.  Finally, 
note  the  convention  used  in  bifurcation  diagrams  is  to  identify  stable 
equilibrium  branches  with  solid  lines  and  unstable  equilibrium  branches 
with  broken  or  dashed  lines;  other  graphical  conventions  will  be  brought 
to  the  reader's  attention  as  required. 
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Figure  5.  Bifurcation  Diagram  for  Example 
Problem 


The  bifurcation  diagram  now  shows  part  of  its  value  in  that  the 
local,  and  at  times  global,  behavior  of  trajectories  can  be  predicted. 
For  this  example  problem  trajectories  in  the  vicinity  (domain  of 
attraction)  of  a  point  on  the  stable  equilibrium  branches  will  oscillate 
about  the  equilibrium  point  since  no  damping  is  included.  If  damping  were 
added  the  trajectories  would  be  attracted  to  the  equilibrium  point.  This 
can  be  visualized  in  thinking  of  the  mass  at  an  angle  9g  that  corresponds 
with  the  stable  branch  for  a  given  value  of  X.  Given  a  small  perturbation 
in  any  direction  the  mass  will  return  to  the  equilibrium  point,  that  is  it 
will  be  attracted  to  the  stable  branch.  Trajectories  in  the  vicinity 
(domain  of  repulsion)  of  a  point  on  one  of  the  unstable  solution  branches 
will  be  repelled.  Again  think  of  the  pendulum  mass  at  an  angle  ©e  that 
corresponds  to  a  point  on  one  of  the  unstable  branches.  Given  a  small 
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perturbation  the  mass  will  not  return  to  its  original  position  or 
trajectory  but  will  move  away  and  in  the  example  problem  will  swing  around 
the  shaft.  Note  finally,  that  the  domain  of  attraction  can  be  limited. 
In  the  example  problem,  even  if  the  mass  is  in  a  stable  position  in  the 
lower  part  of  the  pendulums  arc,  if  given  sufficient  perturbation  in  the 
direction  of  the  force  F  the  mass  can  be  made  to  swing  around  the  shaft 
and  in  the  absence  of  any  damping  action  will  not  settle  back  to  a  stable 
position;  i.e.  the  domain  of  attraction  is  limited.  Finally,  note  for 
values  of  X  >  1  and  A  <  -1  the  pendulum  will  swing  continuously  about  the 
shaft  -  -  no  equilibrium  exists. 

The  foregoing  example  has  provided  a  good  picture  of  the  very  basics 
of  bifurcations  analysis  and  some  of  the  phenomena  that  may  occur  when 
analyzing  nonlinear  problems.  In  addition  it  has  provided  some  physical 
significance  to  limit  points.  It  was  however  limited  in  its  ability  to 
fully  demonstrate  possible  behavior  as  it  has  only  one  degree  of  freedom. 
Problems  with  three  or  more  degrees  of  freedom  can  develop  many  other 
phenomena  and  a  variety  of  ways  to  exchange  stability. 

Limit  Cycles  (Orbits) 

Hopf  bifurcations,  which  are  a  central  feature  in  the  study  at  hand, 
give  rise  to  periodic  solution  branches  which  on  a  bifurcation  diagram  are 
shown  by  plotting  the  maximum  and/or  minimum  values  of  limit  cycles 
(Figure  6).  These  surfaces  projected  in  the  phase  plane,  or  phase  space 
for  systems  with  degrees  of  freedom  greater  than  two,  are  called  limit 
cycles  or  orbits  and  represent  a  surface  of  periodic  solutions  surrounding 
a  equilibrium  solution  branch  (Figure  7) . 
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Figure  6.  Limit  of  Periodic  Motion  about  a 
Solution  Branch 
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Just  as  with  equilibrium  solution  branches,  periodic  solution 
branches  can  be  stable  or  unstable  and  will  either  attract  or  repel 
trajectories  within  their  domain  of  influence  (12:12-25).  Local  stability 
of  the  periodic  solutions  is  based  on  Floquet  theory  involving  the 
monodromy  matrix  (12:240-248).  For  the  purpose  of  this  study  one  need 
only  know  the  relationship  of  the  eigenvalues  of  the  monodromy  matrix  with 
the  criteria  for  local  stability  of  the  periodic  solution  branch.  The 
monodromy  matrix  always  has  one  eigenvalue  (Floquet  multiplier)  in  the 
complex  plane  at  z=l.  This  is  subsequently  used  as  a  test  for  accuracy  in 
calculating  the  other  multipliers  in  AUTO.  Through  establishing  a 
relationship  between  the  remaining  Floquet  multipliers  and  a  Poincard  map 
or  return  map  it  is  determined  that  if  the  modulus  of  the  remaining 
Floquet  multipliers  are  each  less  than  unity  then  the  periodic  solution  is 
stable  (i.e.  attracting).  AUTO  computes  the  Floquet  multipliers  and  uses 
this  to  determine  stability.  As  mentioned  above  the  accuracy  of  this 
calculation  must  be  checked  based  on  the  one  Floquet  multiplier  at  2=1. 
Since  five  ordinary  differential  equations  make  up  the  full  set  of 
equations  needed  to  model  the  longitudinal  dynamics  for  this  study  then 
for  the  periodic  branches  to  be  stable  four  Floquet  multipliers  must  each 
have  a  modulus  less  than  unity  and  one  must  have  a  modulus  equal  to  unity. 
Further  discussion  of  nonlinear  behavior  is  beyond  that  required  to  begin 
the  analysis,  these  concepts  will  be  expanded  on  as  the  results  of  the 
study  are  examined. 
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Two  Parameter  Continuation: 


Upon  locating  a  Hopf  bifurcation  point  or  limit  point  AUTO  can  be 
used  to  perform  a  two  parameter  continuation  which  will  show  the  evolution 
of  bifurcation  points  as  the  second  parameter  is  varied;  ordinary 
equilibria  is  ignored.  For  the  study  here  the  two  parameter  continuation 
will  be  used  to  determine  the  value  of  a  gain  in  a  feedback  loop  to 
attempt  to  stabilize  a  system.  Specifically  by  plotting  one  parameter 
against  another  a  curve  of  limit  points  or  Hopf  points  is  generated.  If 
for  instance,  the  first  parameter  is  the  body  flap  deflection  or  a  pitch 
command  and  the  second  parameter  is  the  pitch  rate  feedback  gain  in  a 
pitch  attitude  feedback  loop,  the  plot  of  the  two  parameters  would  show 
how  the  Hopf  point  (the  point  of  stability  exchange)  moves  with  the  pitch 
rate  gain.  A  typical  plot  might  look  something  like  Figure  8. 

In  this  figure  one  can  see  that  for  values  of  pitch  rate  gain  the  Hopf 
point  will  move  along  the  curve.  If  one  wishes  to  delay  the  occurrence  of 
the  Hopf  point  relative  to  the  magnitude  of  the  body  flap  deflection 
(thereby  delaying  the  exchange  of  stability)  one  would  raise  or  lower  the 
gain  as  appropriate  until  the  Hopf  point  is  moved  far  enough  to  meet  the 
system  requirements  based  on  the  body  flap  deflection.  Some  simple 
feedback  techniques  are  examined  further  in  Section  IV. 
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Figure  8 


Curve  of  Hopf  Points  for  Two  Parameter 
Continuation  with  Pitch  Rate  Gain  and 
Body  Flap  Deflection 
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III.  Bifurcation  Analysis  of  Longitudinal  Dynamics 


To  begin  the  actual  analysis  the  set  of  ordinary  differential  equations 
given  by  eqns(7,9,10,13  and  16)  must  be  solved  for  an  equilibrium  point. 
As  stated  before  this  equilibrium  starting  solution  is  required  for  the 
continuation  method  to  begin.  Taking  the  equations  of  motion  and 
substituting  the  force  and  moment  coefficients  yields: 


T  ,  v  PV2SCD  .  ,  . 

V  =  —  cos (a)  -  — -  -  gsin(y) 

m  2m 


(30) 


4  •  1  -  -~sin(a)  -  ^  +  SM-cosiy) 

mV  2m  i 


(31) 


5  =  g  +  —cos  (y) 
r 


(32) 


=  P^.c.  -  |J[[i£li£]sin(20) 

2  ly  2  *{  Jy  ) 


(33) 


f  =  Vsin(y) 


(7) 


where  the  force  and  moment  coefficients  are  as  defined  before  in  eqns(20, 
21  and  22)  and  are  again  repeated  below  for  convenience. 


+  Cx j,  a2 


(20) 


~  C L0  +  Q*  0 


(21) 
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CjB  =  ^mo  +  ®  +  ^XXJ  + 


(22) 


The  equilibrium  values  are  obtained  by  first  setting  the  left  hand  side  of 
equations  (7,  30,  31,  32  and  33)  to  zero  and  solving  for  the  states  (V,  a, 
q,  6,  and  r) .  To  make  things  easy  y  is  set  to  zero  (y  =  0),  which  means 
a  =  9  =  0  .  From  this  point  the  following  equations  are  obtained  that 
describe  the  requirements  for  equilibrium  with  r  =  rc,  a  =  9  =  0: 


Yo  =  0-0 


(34) 


[  9(r0)r0 

1+ 

fp(r0)ar0C^ 

l  2  m  j 

So  =  - 


Yi 
1 0 


(35) 


(36) 


T(r0 )  =  -|p  {r0)SCDvt 


(37) 


Cm  =  0.0 

in 


(38) 


This  set  of  equations  is  programmed  into  the  user  provided  subroutine  FUNC 
which  AUTO  requires  to  find  the  flow  of  equilibrium  solutions  as  a 
specified  parameter  is  varied.  Note  the  body  flap  parameter  (5bf)  shows 
up  in  the  last  term  of  eqn(22).  Also  of  interest  is  the  value  of  the 
thrust  at  the  starting  equilibrium  point  -  -  it  equals  the  drag  as  was 
discussed  earlier  for  the  thrust  laws.  The  AUTO  software  package  was  used 
to  develop  the  equilibrium  solution  branches  and  identify  bifurcation 
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points  for  the  system  described  by  eqns(7,  30,  31,  32  and  33).  The 
parameters  used  to  conduct  the  analysis  were  a  body  flap  parameter  (6bf) 
and  a  parameter  which  scaled  the  thrust  (6T);  simulating  a  throttle.  This 
study  was  conducted  for  several  starting  equilibrium  altitudes,  which 
amounts  to  selecting  a  starting  altitude  (r  =  rj),  and  setting  the  angle 
of  attack  and  pitch  angle  to  zero  and  letting  AUTO  solve  for  the  velocity, 
pitch  rate  and  thrust.  The  starting  altitudes  ranged  from  50,000  ft  to 
700,000  ft.  Again  note  that  for  each  starting  equilibrium  altitude  a 
different  value  of  the  initial  equilibrium  thrust  is  obtained  since,  as 
discussed  before,  thrust  must  equal  drag  at  the  starting  equilibrium 
point.  Note  that  the  thrust  decreases  with  increasing  altitude.  The 
resulting  bifurcation  diagrams  are  identified  by  their  starting  altitude 
as  the  thrust  level,  which  is  fixed  by  the  drag  at  the  starting  altitude, 
makes  a  difference  in  the  behavior  of  the  system. 

One  can  see  the  equations  of  motion  are  clearly  nonlinear  with  the 
states  all  interrelated,  However  two  states  exert  the  primary  influence 
by  virtue  of  the  type  of  problem  being  analyzed,  these  are  the  velocity 
and  the  radius.  It  is  the  way  these  two  states  vary  to  achieve 
equilibrium  and  the  fact  that  the  behavior  being  observed  is  the  behavior 
of  the  equilibrium  solution  path  that  makes  for  results  that  are  not 
intuitively  obvious  relative  to  the  effects  of  changing  the  body  flap  or 
throttle  (or  more  precisely  thrust  variations).  This  point  must  be 
emphasized  since  most  traditional  dynamists  are  used  to  dealing  with  the 
time  history  of  trajectories  or  frequency  response  given  some  control 
input.  As  discussed  in  section  II  the  bifurcation  diagrams  which  are  used 
here  are  not  time  histories  but  a  collection  of  equilibrium  points  that 
provide  the  value  of  the  states  relative  to  a  parameter. 
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The  fact  that  the  solution  branch  is  made  up  of  equilibrium  points  implies 
that  no  change  of  the  state  variables  occurs  at  the  individual  points 
along  the  curve;  this  can  at  times  cause  confusion  so  beware. 

Three  cases  were  initially  considered  based  on  the  three  separate 
thrust  laws  discussed  in  previous  sections.  The  results  of  the  rocket 
whose  thrust  varied  with  altitude  differed  little  from  the  more  standard 
constant  thrust  rocket  used  by  most  previous  authors  working  on  this  topic 
so  the  case  was  dropped.  The  remaining  two  cases  were  studied  to 
investigate  the  behavior  the  rather  simple  nonlinear  model  would  generate. 
What  follows  is  a  discussion  of  the  results.  Since  the  effect  of  the  body 
flap  and  throttle  are  significantly  different  it  is  best  to  discuss  them 
separately. 

Body  Flap  Parameter  ( 5bf )  Variation 

The  body  flap  parameter  mathematically  represents  a  deflection  of 
the  body  flap  in  degrees.  This  value  is  changed  to  radians  and  affects 
the  value  of  CN  via  eqn{22).  A  change  in  the  value  of  CN  generally  causes 
a  change  in  the  vehicles  angle  of  attack  and  pitch  attitude.  It  is 
interesting  to  note  that  the  primary  influence  of  the  body  flap  on  the 
behavior  of  the  equilibrium  solution  path  is,  that  in  changing  the  angle 
of  attack  the  value  of  the  lift  coefficient  is  changed  which  is  a  key 
parameter  in  establishing  the  value  of  the  velocity  for  a  equilibrium 
orbit  (16:321-344).  Therefore  it  is  best  to  think  of  the  body  flap  as  a 
control  by  which  lift  is  modulated. 

To  understand  the  equilibrium  solutions  obtained  relative  to  the 
body  flap  parameter  one  must  examine  the  relationship  of  the  velocity  and 
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radius  with  C-b.  The  following  equation  relates  the  velocity  to  and 
radius  and  is  central  in  understanding  what  is  occurring  when  looking  at 
a  equilibrium  solution  path  displayed  in  a  bifurcation  diagram. 
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This  equation  results  from  the  homogenous  set  of  equations  for  equilibrium 
and  represents  the  velocity  -  altitude  relationship  for  everything  from  an 
unpowered  satellite  or  lifting  vehicle  in  equilibrium  orbit  to  the  case 
here  of  a  powered  lifting  vehicle  in  equilibrium  orbit.  The  velocity 
plays  a  key  role  in  providing  forces  sufficient  to  balance  the  weight  of 
the  vehicle.  At  high  altitudes  the  centrifugal  force  is  the  primary 
means  by  which  the  vehicle  balances  the  weight,  where  at  lower  altitudes 
the  lift,  which  is  a  function  of  velocity  is  the  primary  balance  to  the 
weight.  The  way  m  which  the  equilibrium  solution  path  moves  in  order  to 
balance  all  forces  is  quite  interesting  and  as  stated  before  not  always 
obvious.  Those  interested  in  knowing  more  about  how  velocity,  altitude 
and  the  lift  coefficient  are  related  in  orbital  flight  are  refered  to 
reference  (16) . 

The  investigation  was  conducted  by  performing  a  continuation  from 
the  equilibrium  starting  point  (i.e.  r=r5)  using  the  body  flap  parameter 
to  control  the  initial  direction  of  the  continuation.  A  body  flap  sweep 
is  defined  as  the  summation  of  the  equilibrium  branches  obtained  from 
performing  the  continuation  in  the  directions  associated  with  a  positive 
flap  deflection  and  a  negative  flap  deflection  (control  surface  movement 
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downward  being  positive) .  The  flap  was  constrained  to  +/-  90  degrees  by 
fixing  the  limiting  values  for  the  parameter  in  AUTO.  Starting  altitudes 
for  the  constant  thrust  rocket  case  were  varied  between  100,000  ft  and 
700,000  ft  and  for  the  air-breathing  engine  case  ranged  from  50,000  ft  to 
400,000  ft.  Note  in  all  cases  as  the  starting  altitude  is  increased  the 
starting  equilibrium  value  of  the  thrust  is  decreased.  The  equilibrium 
solution  paths  were  tabulated,  stability  determined  and  all  simple 
bifurcations,  limit  points  and  Hopf  bifurcation  points  were  located. 
Having  mapped  out  the  equilibrium  solution  and  identified  the  various 
singular  or  bifurcation  points  any  Hopf  bifurcations  were  continued  to 
obtain  the  limit  cycles.  This  process  is  accomplished  by  taking  the 
equilibrium  conditions  at  the  point  identified  as  a  Hopf  bifurcation  as  a 
starting  point  for  AUTO's  continuation  method.  Specific  software  routines 

t 

in  AUTO  are  used  to  perform  the  required  functions  to  obtain  the  periodic 
solution  branch.  These  data  are  generally  interpreted  graphically  to 
obtain  a  general  feel  for  the  local  behavior  of  trajectories  in  the 
vicinity  of  the  solution  branches;  these  graphs  are  known  as  bifurcation 
diagrams.  Since  bifurcation  diagrams  are  meant  to  convey  information 
about  the  behavior  of  the  system,  it  seems  only  natural  that  a  method  or 
convention  be  established  for  presenting  data  on  these.  The  reader  is 
encouraged  to  take  note  of  the  following  rules  for  conveying  information 
about  the  types  of  solution  branches  and  their  local  stability 
characteristics. 

Equilibrium  solution  branches  are  presented  as  lines.  Solid  lines 
indicate  stable  solution  branches  and  any  type  of  broken  or  dashed 
lines  indicate  an  unstable  equilibrium  solution  branch. 

Periodic  solution  branches  are  shown  as  circles  or  dots  which 
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generally  indicate  the  maximum  amplitude  of  the  periodic  motion. 
Stable  periodic  solution  branches  are  shown  as  solid  dots  or  filled 
circles.  Unstable  periodic  solution  branches  are  indicated  by  open 
circles. 

These  conventions  are  adhered  to  throughout  this  paper. 

Constant  Thrust  Rocket  Case.  The  following  bifurcation  diagrams 

were  generated  using  AUTO  as  described  previously,  however  further 

explanation  of  the  way  the  constant  thrust  value  is  determined.  As  stated 

several  times  before,  and  shown  explicitly  in  eqn(37),  the  thrust  equals 

the  drag  at  the  starting  equilibrium  solution  with  y  =  0  radians.  For 

the  constant  thrust  rocket  case  this  value  obviously  is  fixed  over  the 

entire  body  flap  sweep.  This  means  that  two  equilibrium  solution  points 

with  the  same  altitude  but  obtained  from  body  flap  sweeps  starting  at  two 

different  altitudes  will  not  have  the  same  value  of  thrust  and  in  general 

will  have  different  values  for  the  other  states  as  well.  Note  finally,  as 

starting  altitudes  are  increased  the  value  of  the  thrust  decreases. 

Figures  9  and  10  show  a  collection  of  bifurcation  diagrams  for  each 

state  (note  a  =  0  )  for  the  body  flap  (5„<)  sweeps  from  100,000  and  300,000 

•»  * 

ft  respectively.  While  the  behavior  is  nonlinear  there  is  not  a  great 
deal  of  interest  occurring  over  most  of  the  equilibrium  branches. 

Figure  9  is  characteristic  of  the  constant  thrust  rocket  case  with  5^  as 
the  parameter  for  starting  altitudes  less  than  150,000  ft  and  Figure  10  is 
characteristic  of  the  constant  thrust  rocket  case  with  5^  as  the  parameter 
for  starting  altitudes  between  150,000  ft  and  about  360,000  ft.  One  of 
the  first  things  to  note  in  each  figure  is  the  system  is  unstable 
(indicated  by  the  dashed  line). 
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Figure  9.  Bifurcation  Diagrams  for  Body  Flap  Sweep  from  100  kft 
Constant  Thrust  Rocket 


-50  0  50  -50  0  50 

Body  Flap  Deflection  [deg]  Body  Flap  Deflection  [deg] 


Body  Flap  DeHection  [deg]  Body  Flap  Deflection  [deg] 

Figure  10.  Bifurcation  Diagrams  for  Body  Flap  Sweep  from  300  kft 
Constant  Thrust  Rocket 


The  instability  is  caused  by  the  nonoscillatory  height  mode  and  was  the 
case  for  both  thrust  laws  representing  a  rocket.  This  behavior  is 
consistent  with  the  previous  studies  and  follows  as  a  result  of  the  way 
the  sum  of  the  longitudinal  forces  change  with  respect  to  altitude 
(13:472;  7:283).  Interestingly  enough  the  height  mode  would  generally 
stabilize  at  some  point  for  starting  altitudes  between  about  360,000  ft 
and  530,000  ft  and  would  occur  associated  with  a  limit  point. 
Frustratingly  this  is  about  the  point  where  the  mode  normally  thought  of 
as  the  phugoid  mode  (based  on  the  longer  period)  would  go  unstable. 

From  the  stand  point  of  nonlinear  analysis  not  much  of  interest  is 
occurring  for  the  body  flap  sweeps  with  starting  altitudes  less  than  about 
360,000  ft.  Specifically,  the  system  is  unstable  with  generally  no  limit 
points  from  which  jump  phenomena  may  occur  and  no  simple  bifurcations. 
Only  if  the  continuation  of  a  equilibrium  branch  associated  with 
increasing  altitude  is  allowed  to  go  long  enough,  to  where  the  aerodynamic 
pitch  damping  is  lost  due  to  the  very  low  density  at  very  high  altitudes, 
will  a  Hopf  point  be  found.  For  the  body  flap  sweep  from  300,000  ft 
(Figure  10)  the  Hopf  point  occurs  at  5^  =  +/-  57.9°  ,  and  an  altitude  of 
615,120  ft  for  the  branch  associated  with  a  negative  body  flap  deflection 
and  614,970  ft  for  the  other  branch;  this  symmetric  behavior  is  not  as 
closely  followed  at  lower  starting  altitudes  where  the  thrust  levels  and 
densities  are  greater. 

Concentrating  now  on  the  unstable  periodic  branches  found  from  the 
body  flap  sweep  from  300,000  ft  one  can  see  from  the  bifurcation  diagrams 
in  Figure  10  that  as  the  unstable  periodic  branch  progresses  the  amplitude 
of  the  limit  cycles  become  fairly  substantial,  however  the  period  is  on 
the  order  of  5000  seconds  and  since  the  periodic  branch  is  unstable 
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trajectories  would  not  approach  it.  Recall  that  the  limit  cycles  of  a 
periodic  branch,  indicated  by  circles,  show  the  maximum  and/or  minimum 
amplitude  of  the  ■;  *  '.die  motion.  The  stability  of  the  periodic  branches 
is  interesting  in  that  it  is  just  barely  unstable  with  a  conjugate  pair  of 
Floquet  multipliers  just  outside  the  unit  circle.  The  height  mode  at  the 
point  where  the  phugoid  mode  goes  unstable  is  very  near  the  imaginary  axis 
with  a  time  constant  on  the  order  of  10^  seconds;  this  is  consistent  with 
previous  studies  (4:786). 

Figure  11  shows  the  time  history  over  one  period  for  the  right 
periodic  branch  with  5^-  =  56.73°  and  5^=  42.25°.  Figure  12  shows  the 
time  history  over  one  period  for  the  left  periodic  branch  with  5hf  =  -56.6° 
and  6„;  =  -41.8°.  Note  that  near  the  Hopf  point  on  either  branch  the 
motion  is  slight  (ie  just  leaving  equilibrium)  but  as  the  parameter  is 
changed  to  move  along  the  periodic  branch  the  motion  increases.  Once 
again  the  behavior  of  the  periodic  branches  shown  in  this  body  flap  sweep 
(from  300,000  ft)  is  characteristic  of  the  behavior  of  the  periodic 
branches  occurring  from  body  flap  sweeps  starting  at  "lower"  altitudes 
that  subsequently  extend  to  altitudes  above  500,000  ft  where  pitch 
stability  is  lost. 

In  order  to  provide  a  complete  look  at  the  periodic  behavior  of  the 
limit  cycles,  as  well  as  provide  a  connection  with  more  classic 
longitudinal  analysis,  the  limit  cycles  are  projected  into  the  phase  plane 
with  the  flight  path  angle.  Figure  13  shows  two  limit  cycles  from  the 
right  periodic  branch;  one  for  6^  =  56.7°  and  one  for  5^  =  42.25°. 

Figure  14  shows  two  limit  cycles  from  the  left  periodic  branch;  one  at 
6.,-  =  -56.64°  and  one  for  6^  =  -41.8°.  In  looking  at  these  figures  one 
sees  clearly  that  motion  is  associated  with  the  periodic  branches. 
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Figure  11.  Left  Branch  Limit  Cycles  for  6^  =  -56.6°and  -41.8° 
Body  Flap  Sweep  from  300  kft  :  Const  Thrust  Rocket 


Figure  12.  Right  Branch  Limit  Cycles  for  6^  =  56.73°  and  42.25° 
Body  Flap  Sweep  for  300  Kft  :  Const  Thrust  Rocket 
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Figure  13. 


Right  Branch  Limit  Cycles  in  Phase  Plane  for  =  56.73° 
and  42.25°.  Sweep  from  300  kft  :  Const  Thrust* Rocket 
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Figure  14.  Left  Branch  Limit  Cycles  in  Phase  Plane  for  5^  =  -56.64° 
and  -41.8°.  Sweep  from  300  kft  :  Const  Thrust  Rocket 
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Also,  the  behavior  of  the  translational  states  is  very  much  like  the 
classic  phugoid  mode,  i.e.  showing  a  steady  exchange  of  potential  and 
kinetic  energy.  While  the  periodic  branches  are  unstable,  the  growth  of 
the  nonlinear  behavior  in  velocity  and  altitude  remains  sinusoidal. In 
terms  of  the  rotational  states,  one  sees  somewhat  more  complex  behavior 
which  should  be  expected  since  the  periodic  branches  arose  from  the  loss 
of  pitch  stability  at  the  very  high  altitude. 

The  most  interesting  body  flap  sweep  for  the  constant  thrust  rocket 
case  resulted  from  the  starting  altitude  of  400,000  ft.  The  body  flap 
sweep  bifurcation  diagrams  and  expanded  views  for  a  and  altitude  are  shown 
in  Figures  15,  16  and  17.  The  behavior  of  the  limit  cycles  of  the 

periodic  solution  branches  are  certainly  visually  interesting.  Note  the 
periodic  branches  contain  several  limit  points  which  explains  their 
complex  twisting  about. 

The  Hopf  bifurcation  occurred  at  5b£  "  +/-  8.93°  and  generally 

speaking  is  not  of  great  significance  since  the  equilibrium  branch  was 
unstable  to  begin  with  and  the  periodic  branch  starts  out  unstable  and 
encircles  an  unstable  equilibrium  branch.  On  closer  inspection  of  the 
periodic  branches  one  will  see  (Figures  16  and  17)  that  there  is  a  portion 
of  both  periodic  branches,  starting  at  =  +/-  7.82°  and  continuing  to 

5bf  =  +/,“  6-05°'  that  gains  stability  by  the  crossing  of  a  conjugate  pair 

of  Floquet  multipliers.  This  type  of  stability  exchange  is  associated 
with  bifurcation  to  a  torus.  What  this  implies  is  that  trajectories 
within  the  domain  of  attraction  of  the  periodic  branch  will  be  drawn  into 
periodic  motion  with  two  frequencies;  one  describing  the  component  of  the 
motion  in  the  circumference  direction  of  the  torus  and  one  around  the 
cross  section  of  the  torus.  The  path  of  the  trajectory  can  be  visualized 
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Figure  15.  Bifurcation  Diagrams  for  Body  Flap  Sweep  from  400  kft 
Constant  Thrust  Rocket 


as  spiralling  around  the  inside  of  an  inner  tube  (i.e.  torus)  in 
hyperspace  (12:263,264). 

The  accuracy  of  the  solution  is  very  good  with  the  Floquet 
multiplier  that  is  supposed  to  be  equal  to  one  (z=l) ,  precisely  equal  to 
one.  The  fact  that  stability  is  gained  then  lost  so  quickly  indicates 
that  one  or  more  Floquet  multiplier (s)  is (are)  very  near  the  edge  of  the 
unit  circle  and  upon  inspection  of  the  output  from  AUTO  one  finds  this  to 
be  the  case  with  one  pair  inside  the  unit  with  modulus  =  0.99884  and  a 
second  pair  with  modulus  =  0.98.  This  would  indicate  a  weakly  attracting 
limit  cycle  for  the  range  of  5^  where  the  periodic  branches  are  stable. 
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6bf  Sweep  from  400  Kft  :  Const  Thrust  Rocket 


5M'--Sweep  from  400  Kft  :  Const  Thrust  Rocket 


Figure  17.  h  Bifurcation  Diagram  for  Body  Flap  Sweep  from  400  kft 
Constant  Thrust  Rocket 


Figure  18  shows  the  motion  of  the  limit  cycle  over  one  period  for 
periodic  branch  1  (left  branch)  for  a  point  in  the  stable  region  at 
Sjjj  =  -6.83°  and  a  point  near  the  end  of  the  calculated  portion  of  this 
periodic  branch  at  5^  =  -2.89°.  Figure  19  shows  the  motion  of  the  same 
two  limit  cycles  (5^  =  -6.83°  and  5bf  -  -2  .89°)  over  one  period  in  the 
phase  plane  with  the  flight  path  angle.  Figure  20  shows  the  motion  of  the 
limit  cycle  over  one  period  for  periodic  branch  2  (right  branch)  with 
5b£  =  8.73°  in  the  stable  region  and  for  5bb  =  3.57°  near  the  end  of  the 
calculated  branch.  Figure  21  shows  the  limit  cycles  at  6bb  =  8.73°  and 
5^  =  3.57°  in  the  phase  plane  with  flight  path  angle.  Observing  the 
behavior  of  the  limit  cycle  over  one  period  at  several  points  like  this 
shows  why  the  nonlinear  behavior  associated  with  period  solutions  are  so 
interesting.  Looking  at  Figures  18  and  20  one  can  see  a  kind  of  wave 
changing  in  amplitude  as  the  parameter  is  varied.  Figure  22  shows 
qualitatively  the  growth  of  the  nonlinear  behavior  in  a  as  periodic 
branch  2  grows.  It  is  this  type  of  behavior,  for  systems  with  three  or 
more  degrees  of  freedom,  that  leads  to  more  fascinating  subjects  like 
chaos  and  the  Hopf  point  is  as  Seydel  puts  it,  "the  door  which  opens  from 
the  small  room  of  equilibria  to  the  large  hall  of  periodic  solutions 
(12:61)." 

On  a  somewhat  different  note,  the  altitude  where  the  Hopf  point 
occurs  or.  both  branches  is  just  about  450,000  ft.  This  is  in  the  altitude 
range  where  Etkin  determined,  and  later  others  modified,  that  the  period 
of  phugoid  and  pitching  modes  came  very  close  to  each  other  (4:787-788; 
15:17-20).  The  general  conclusion  of  these  earlier  works  is  that  there 
would  be  significant  coupling  of  the  two  modes  at  this  so  called 
"resonance  altitude"  (15 ; 7) . 
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18.  Limit  Cycles  for  Periodic  Branch  1  (5^  =-6.826°  &  -2.89°) 
Body  Flap  Sweep  from  400  kft  :  Const  Thrust  Rocket 


.575 
2.57 
2.565 
2.56 
2.555 

-0.2  -0.1  0  0.1  0.2 
flight  Poth  Angl#  [d#g] 

0.05 
0 

•0.05 
-0.1 
■0. 


-0.2  -0.1  0  0.1  0.2  -0.2  -0.1  0  0.1  0.2 
night  Poth  Angl#  [d#gj  night  Poth  Angl#  (d#g] 

Figure  19.  Limit  Cycles  in  Phase  Plane,  Periodic  Branch  1 

<6l£  =  -6.826°  &  -2.89°)  :  Constant  Thrust  Rocket 
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Figure  20.  Limit  Cycles,  Periodic  Branch  2,(5^  =  8.738°  &  3.57°) 
Body  Flap  Sweep  from  400  kft  :  Const  Thrust  Rocket 


Figure  21.  Limit  Cycle  in  Phase  Plane,  Periodic  Branch  2  (5^  =  8.738° 
&  3.57°),  Sweep  from  400  kft  :  Const  Thrust  Rocxet 
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Figure  22.  Growth  of  Nonlinear  Waveform,  Periodic  Branch  2, 

a  Bifurcation  Diagram,  Body  Flap  Sweep  from  400  kft 


In  looking  at  the  limit  cycles  in  the  preceding  figures  there  does 
not  seem  to  be  the  strong  coupling  predicted.  Some  coupled  motion  is 
evident  in  that  the  rotational  states  (a,  0,  and  q)  go  through  sub¬ 
oscillations  in  each  overall  period  while  the  translational  states 
maintain  a  sinusoidal  motion  with  a  very  regular  period.  Since  the  period 
is  nearly  that  of  the  circular  orbit  for  the  same  geocentric  radius,  it 
would  seem  that  what  is  observed  is  a  barely  unstable  elliptical  orbit 
with  the  vehicle  pitching  about  its  y-axis  at  some  rub-frequency  greater 
than  the  frequency  of  the  orbit  (i.e.  overall  frequency  of  the  limit  cycle 
for  the  given  parameter). 

As  a  final  note,  notice  the  limit  points  on  the  equilibrium  branch 
where  the  equilibrium  solution  path  changes  direction. 


It  is  interesting  in  that  it  exists  and  is  associated  with  the  height  mode 
stabilizing.  Notice  how  the  direction  of  the  body  flap  deflection  changes 
in  order  to  maintain  the  original  direction  of  the  equilibrium  solution 
path. 


Air-breathing  Engine  Case.  The  procedure  for  continuation  and 
subsequent  analysis  for  the  air-breathing  case  with  the  body  flap  as  the 
parameter,  was  the  same  as  that  for  the  constant  thrust  rocket  analysis. 
The  results  are  somewhat  more  interesting  in  that  the  height  mode  is 
generally  stable  for  the  air-breathing  case  below  approximately 
380,000  ft  thus  the  entire  system  is  stable  at  these  "lower"  altitudes. 

A  phenomenon  of  little  physical  significance,  but  interesting 
nonetheless  for  the  air-breathing  case  with  body  flap  sweeps  starting  from 
equilibrium  points  between  approximately  100,000  ft  up  to  approximately 
360,000  ft  is  that  the  velocity  goes  very  nearly  to  zero  for  the  portion 
of  the  body  flap  sweep  that  has  a  negative  body  flap  deflection.  Before 
the  velocity  actually  gets  to  zero,  a  Hopf  bifurcation  occurs,  then  a 
simple  bifurcation  which  has  two  solution  branches.  Of  the  two  branches 
one  stable  and  back-track  the  original  equilibrium  solution  branch  for  all 
the  states,  and  the  other  is  unstable.  The  unstable  branching  solution 
back-tracks  a,  and  altitude,  but  takes  the  negative  of  it's  original  value 
for  velocity  and  pitch  rate.  Figure  23  shows  this  bifurcation  diagram  for 
the  equilibrium  branches  only.  Further  work  is  needed  to  finish  exploring 
this  behavior.  Notice  that  in  Figure  23  the  +  -  symbol  indicates  the 
stable  branching  solution  and  the  x  -  symbol  indicates  the  unstable 
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branching  solution.  Finally  the  original  branch  also  turns  back  on  itself 
after  going  to  nearly  zero  subsequent  to  the  Hopf  point;  this  back¬ 
tracking  original  branch  turns  back  as  an  unstable  branch  but  regains 
stability  as  it  repasses  the  Hopf  point  (as  expected) . 


Figure  23.  Bifurcation  Diagram  for  Body  Flap  Sweep  from  100  kft 
Air-Breathing  Engine 


For  ease  of  discussion,  Figure  24  shows  all  of  the  positive  body  flap 
sweep  from  100,000  ft,  but  only  the  stable  portion  of  the  negative  body 
flap  sweep  from  100,000  ft  for  the  air-breathing  engine  case.  This 
diagram  contains  the  Hopf  bifurcation  but  none  of  the  "back-tracking" 
solutions.  A  big  difference  from  the  previous  case  is  readily  apparent, 
the  system  is  stable  up  to  the  Hopf  point  at  which  time  the  phugoid  mode 
loses  stability  and  note  how  low  the  altitude  is  (  approximately  73,000 
ft).  Notice  this  is  a  subcritical  bifurcation  (12:72);  that  is  an 
unstable  periodic  branch  encircles  a  stable  equilibrium  branch.  What  this 
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implies  is  trajectories  near  the  periodic  solution  branch  but  within  the 
domain  of  attraction  of  the  equilibrium  branch  will  be  drawn  to  the 
stable  equilibrium  branch  and  generally  no  further  changes  of  the  states 
will  occur. 


Figure  24.  Partial  Body  Flap  Sweep  from  100  kft  for  the  Air-Breathing 
engine  case 


Figure  25  shows  the  limit  cycles  over  one  period  for  a  point  near 
the  Hopf  bifurcation  {5^  =  -52.77°  ;  T  =  77  sec)  and  the  point  on  the  end 
of  the  calculated  periodic  branch  (5b{  =  -23.54°  ;  T  =  140  sec).  Clearly 
visible  is  the  increase  in  nonlinear  behavior  as  one  moves  along  the 
periodic  branch.  Near  the  Hopf  point  one  can  see  the  motion  is  nearly 
constant  and  what  little  variation  is  there  is  sinusoidal.  For  points 
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Figure  25.  Limit  Cycles  for  =  -52.77°  &  6^  =  -20.54°,  Body  Flap 
Sweep  from  100  kft  :  Air-Breathing  Engine 


farther  along  the  periodic  branch  the  nonlinear  behavior,  due  to  the 
underlying  nonlinear  equations,  truly  begins  to  blossom. 

Figure  26  shows  the  limit  cycles  for  5^  =  -52.77°  and  -23.54°  in  the 
phase  plane.  Once  again  the  classic  phugoid-like  behavior  of  the 
translational  states  is  seen.  Note  as  before  in  the  constant  thrust  case 
the  rotational  states  have  this  sub-oscillation.  In  contrast  to  what  was 
seen  in  the  constant  thrust  rocket  case,  here  the  phugoid  mode  has  gone 
unstable  at  relatively  low  altitude  (approximately  73,000  ft). 

A  body  flap  sweep  from  400,000  ft  for  the  air-breathing  engine  case 
is  shown  in  Figure  27.  The  general  characteristics  for  this  sweep  are  the 
same  as  for  the  constant  thrust  rocket.  However  here  the  periodic  branch 
remains  unstable.  As  with  the  constant  thrust  rocket  case,  each 
equilibrium  branch  has  a  limit  point  where  the  height  mode  stabilizes. 
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night  Poth  Anglo  [dag]  night  Path  Angla  [dag] 

Figure  26.  Limit  Cycles  in  Phase  Plane  (5^  =  -52.77°  &  -23.54°) 
Body  Flap  Sweep  from  100  kft  :  Air-Breathing  Engine 


There  are  also  several  limit  points  found  along  each  periodic  branch;  just 
as  with  the  constant  thrust  rocket  body  flap  sweep  from  400,000  ft. 

Figures  28  and  29  show  the  limit  cycles  over  one  period  and  the 
limit  cycles  in  the  phase  plane  for  the  right  periodic  branch  with 
5fc£  =  -0.00078°  and  0.0118°.  Figures  30  and  31  show  the  limit  cycles  over 
one  period  and  the  limit  cycles  in  the  phase  plane  for  the  left  periodic 
branch  with  =  -0.0046°  and  -0.011°.  These  points  are  as  before,  used 
to  show  the  behavior  of  the  limit  cycle  near  the  Hopf  point  versus  the 
behavior  farther  down  the  periodic  branch.  The  same  basic  behavior  is 
seen  for  the  translational  states  as  discussed  for  previous  cases. 
However  of  interest  is  the  phase  shift  occurring  along  the  left  periodic 
branch.  It  is  interesting  to  see  how  the  3°ft  branch  differs  markedly 
from  the  right  branch  even  though  the  two  look  relatively  symmetric. 
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Figure  27.  Collection  of  Bifurcation  Diagram  for  Body  Flap  Sweep  from 
400  kft  :  Air-Breathing  Engine  Case 


r 


I 


Figure  28.  Limit  Cycles  over  one  Period  5bf=  -0.0007°  and  0.0118° 

Body  Flap  sweep  from  400  Xft  :  Air-Breathing  Engine  Case 


Figure  29.  Limit  Cycles  in  Phase  Plane  for  6bf=  -0.0007°and  0.0118° 

Body  Flap  Sweep  from  400  kf t  :  Air-Breathing  Engine 
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Figure  30.  Limit  Cycles,  Left  Branch, (  6^  =  -0.0046°&  -O.Olf) 
Body  Flap  Sweep  from  400  kft  :  Air-Breathing  Case 


Figure  31.  Limit  Cycles  in  Phase  Plane,  Left  Branch  (5^  =  -0.0046° 
&  -0.011°),  Sweep  from  400  kft  :  Air-Breathing  Engine 
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Throttle  Parameter  (5T)  Variation 


The  throttle  parameter  acts  as  a  multiplier  to  scale  the  equilibrium 
value  of  thrust.  In  the  case  of  the  constant  thrust  rocket  it  scales  the 
value  of  the  thrust  at  the  starting  altitude  which  was  set  equal  to  the 
drag.  For  the  air-breathing  engine  case  the  throttle  parameter  scales  the 
thrust  as  with  the  rocket  case  (i.e.  which  was  set  equal  to  drag  at  the 
starting  altitude) ,  however  for  the  air-breathing  case  the  thrust  varies 
as  a  function  of  altitude  via  T  =  T0  [p(r)/pjj].  A  change  in  thrust  causes 
a  corresponding  change  in  drag  to  maintain  equilibrium  so  once  again  the 
velocity  and  radius  begin  to  play  an  important  role.  However  an 
interesting  characteristic  of  the  equilibrium  solution  paths  where  the 
throttle  is  the  parameter  is  the  lack  of  any  modulation  in  the  coefficient 
of  lift  or  drag;  which  remain  effectively  constant  at  all  throttle 
settings  and  radii. 

The  analysis  using  the  throttle  as  a  parameter  was  done  in  the  very 
same  manner  as  the  body  flap  parameter.  As  before  a  starting  altitude  was 
selected  (with  a  =  9  =  0  radians)  from  which  the  equilibrium  values  for 
the  remaining  states  were  calculated.  From  this  point  the  throttle 
parameter  was  first  increased  from  6T=1.0  then  decreased.  This  yielded 
a  complete  throttle  sweep. 

Constant  Thrust  Rocket  Case.  Not  much  happened  with  this  case.  From 
Figures  32  and  33  one  can  see  the  system  is  unstable,  due  to  the  height 
mode,  and  is  nonlinear  but  no  bifurcations  were  detected.  It  appears 
without  the  pitching  associated  with  the  body  flap  there  is  little  to 
drive  the  nonlinear  nature  of  the  problem.  In  contrast  to  the  cases  where 
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the  body  flap  was  used  as  the  parameter  the  angle  of  attack  in 
cases  is  basically  zero. 


these 


Figure  32.  Collection  of  Bifurcation  Diagrams  for  Throttle  Sweep  from 
100  kft  :  Constant  Thrust  Rocket 
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Figure  33.  Collection  of  Bifurcation  Diagrams  for  Throttle  Sweep  from 
400  kft  :  Constant  Thrust  Rocket 
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Air-Breathing  Engine  Case.  For  the  air-breathing  case  the  most 
interesting  behavior  was  at  altitudes  below  400,000  ft.  The  following 
discussion  characterizes  the  general  nature  of  what  was  found  at  these 
"lower"  altitudes.  Figure  34  shows  the  5T  sweep  for  the  air-breathing 
engine  from  100,000  ft.  For  the  portion  of  the  sweep  where  5T  <  1.0  the 
system  remains  stable.  For  the  portion  of  the  sweep  where  5T  >  1.0  the 
system  loses  stability  subcritically  at  a  limit  point  (5T  =  18.99  )  where 
the  height  mode  crosses  the  imaginary  axis.  Just  after  this  occurs  (5T  = 
18.93)  a  Hopf  Bifurcation  point  is  found.  Since  the  limit  point  preceded 
the  Hopf  point  this  bifurcation  is  not  classified  either  supercritical  or 
subcritical  and  as  before  is  really  of  little  physical  value  other  than  to 
perhaps  give  an  idea  of  the  bound  on  the  allowable  perturbation  to  remain 
in  the  vicinity  of  the  equilibrium  solution  branch. 

Figure  35  shows  an  expanded  view  of  the  area  around  the  Hopf  point 
and  Figure  36  shows  just  the  maximum  limit  cycles  from  the  bifurcation 
diagram.  These  limit  cycles  show  quite  a  variation  of  amplitude  for  the 
rotational  states  along  the  periodic  branch  as  6T  is  increased.  Looking 
at  Figure  37  one  sees  smooth  sinusoidal  behavior  in  the  translational 
states  even  though  the  limit  cycles  are  for  points  well  toward  the  end  of 
the  calculated  portion  of  the  periodic  branch  (5T=  19.063  and  19.078). 
The  rotational  states  however,  display  some  relatively  high  frequency  sub¬ 
oscillations.  Although  the  amplitudes  of  these  sub-oscillations  are  quite 
small  and  the  period  is  on  the  order  of  150  seconds  (see  Figure  38).  In 
examining  the  phase  plane  representations  of  the  limit  cycles  versus  the 
flight  path  angle  (Figures  39  and  40)  one  sees  the  translational  states 
clearly  displaying  the  motion  that  can  be  associated  with  an  elliptical 
orbit.  Further  support  of  this  view  is  given  from  Figure  41,  where  the 


63 


variation  in  the  overall  parameter  dependent  period  of  the  periodic  branch 
is  shown  relative  to  the  circular  orbital  period  for  the  values  of  the 
states  at  the  given  5T. 
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Figure  34.  Collection  of  Bifurcation  Diagrams  for  Throttle  Sweep  from 
100  kft  for  Air-Breathing  Engine  Case 
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Figure  35.  Expanded  View  of  Region  near  Hopf  Point  for  5T  sweep 
from  100  kft  for  Air-Breathing  Engine  Case 
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Figure  36.  Maximum  Limit  Cycles  for  the  5T  sweep  from  100  kft  for 
Air-Breathing  Engine  Case 
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38.  Expanded  View  of  a  Limit  Cycle,  5T  =  19.063  Throttle  Sweep 
from  100  kft  :  Air-Breathing  Engine 
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Figure  39.  Limit  Cycles  in  Phase  Plane  for  5T  =  19.063 
Air-Breathing  Engine  Case  from  100  Kft 
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Figure  40.  Limit  Cycles  in  Phase  Plane  for  5T  =  19.078 
Air-Breathing  Engine  Case  from  100  kft 


Figure  41.  Comparison  of  Circular  Orbital  Period  and  Period  of  Limit 
Cycles  for  Throttle  Sweep  from  100  kft  :  Air-Breathing 
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IV.  Model  Stabilization 


Simple  Feedback  Options 

In  Stengel's  work  (13)  he  presented  a  summary  of  several  possible 
feedback  schemes  using  several  controls.  Table  1  in  section  I  is  a 
reproduction  of  the  table  found  in  Stengel's  paper.  In  looking  at  these 
possible  schemes  one  becomes  aware  of  the  dichotomy  regarding  the  height 
mode  and  the  phugoid  mode.  Most  stabilizing  feedback  for  the  height  mode 
destabilizes  the  phugoid  and  vice  versa.  In  Berry’s  work  (2)  similar 
feedback  options  were  tried  for  the  linear  approximation  technique  and 
found  to  display  the  same  behavior  as  found  by  Stengel.  This  "inverse" 
relationship  between  the  height  mode  and  the  phugoid  mode  plus  the 
restrictions  imposed  by  the  simplicity  of  the  vehicle  model  made  it  beyond 
the  scope  of  this  thesis  to  actually  stabilize  the  height  mode.  The 
success  in  stabilizing  the  height  mode  and  not  destabilizing  the  phugoid 
lay  in  developing  a  control  law/technique  to  properly  modify  the  way  the 
longitudinal  forces  vary  with  height  and  velocity  (as  discussed  in 
sec. ion  I).  Minor  success  was  experienced  m  dealing  with  the  phugoid  by 
using  pitch  rate  feedback  to  the  body  flap.  It  should  be  noted  that  pitch 
rate  feedback  in  general  has  little  afuc-  on  the  phugoid  roots  however  it 
was  a  technique  that  could  be  easily  managed  and  did  demonstrate  the 
concept.  For  a  given  change  in  the  pitch  rate  feedback  the  phugoid  mode 
could  be  improved,  but  so  slightly  that  in  a  practical  sense  it  was 
worthless.  Figure  42  shows  the  basic  feedback  loop  with  the  pitch  rate 
relative  to  the  earth  fed  back  in  a  negative  feedback  loop  to  the  body 
flap  (5b{). 
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Figure  42.  Pitch  Rate  Feedback  Loop  for  Model  Stabilization 


Note  the  parameter  now  becomes  the  value  input  as  the  command  (5^) . 

The  new  value  of  5^  for  use  in  eqn(22),  the  moment  coefficient  equation 
is  given  by  the  following: 

(« -  d)  ^  <»> 

where:  5^  =  new  parameter  for  body  flap  sweep  control  (deg) 

Kq  *  pitch  rate  feedback  gain 

Figure  43  is  sufficiently  representative  of  all  the  cases  where  pitch  rate 
feedback  was  used.  It  must  be  emphasized  that  the  curve  in  Figure  41  is 
a  curve  of  Hopf  bifurcation  points.  What  can  be  seen  is  that  as  the  value 
of  Kq  is  changed  the  location  of  the  Hopf  point  relative  to  the  body  flap 
deflection  is  changed.  In  this  case  a  gain  of  Kq=20  pushes  the  Hopf  point 
down  the  curve  the  farthest.  However,  as  one  can  clearly  see  the 
improvement  is  extremely  small;  to  the  point  of  being  basically 
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constant  thrust  rocket  system  and  attempting  to  use  alternative  feedback 
of  other  states,  it  seems  that  without  recasting  the  mathematical  model  to 
allow  for  more  reasonable  feedback  control,  say  with  attitude  command 
inputs,  significant  stabilizing  routines  cannot  be  obtained. 


Figure  43.  Movement  of  Hopf  Bifurcation  given  Pitch  Rate  Feedback  to 
the  body  flap.  Air-Breathing  Engine  from  100  kft 
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V.  Conclusions  and  Recommendations 


Conclusions 


Overall  this  study  has  not  produced  results  that  would  be  considered 
"Earth  shaking"  or  significantly  different  from  previous  work.  It  has 
however  resulted  in  some  worthwhile  accomplishments,  not  the  least  of 
which  is  demonstrating  how  easy  it  is  to  use  bifurcation  analysis  on 
problems  related  to  hypersonic  flight.  This  method  provided  much  the  same 
information  as  obtained  by  other  authors  using  perturbation  techniques, 
yet  gave  a  much  greater  view  of  the  actual  effects  of  the  nonlinear  nature 
of  the  problem. 

In  terms  of  comparison  to  previous  work,  it  was  found  that  the 
period  and  damping  of  the  phugoid  and  pitching  modes  was  similar  to  the 
behavior  discovered  by  Vihn  and  Dobrzelecki  and  verified  by  Markopoulos, 
et  al.  (15:16-18;  7:286,287).  Their  study  showed  that  the  two  modes  do 
not  cross  for  the  linearized  model  as  Etkin  had  concluded,  but  instead 
come  very  close  together  then  diverge  (4:785,786);  this  was  the  case  here 
as  well.  From  past  work  the  behavior  associated  with  the  two  modes  in  the 
vicinity  of  the  "resonance  altitude"  ( 1 5 ; 1 6 )  are  expected  to  be  coupled 
and  behave  nothing  like  that  expected  cf  the  classic  pitching  and  phugoid 
modes.  It  was  shown  here  that  there  is  significant  departure  from  the 
classic  behavior  of  the  rotational  states,  in  that  they  show  significant 
motion  when  the  phugoid  mode  goes  unstable.  However  the  translational 
states  act  basically  as  expected.  Looking  at  the  limit  cycles  of  the 
periodic  branch  associated  with  the  "resonance  altitude"  one  sees  what 
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could  be  described  as  a  one  way  coupling  from  the  translational  states  to 
the  rotational  states,  with  the  rotational  states  experiencing 
sub-oscillating  for  each  period  of  the  translational  states.  In  addition 
the  loss  of  stability  at  high  altitude  which  Etkin  concluded  would  occur 
due  to  the  loss  of  pitch  damping  and  a  destabilizing  moment  due  to  the 
effect  of  the  gravity  gradient  was  seen  (4:785,787). 

For  the  body  flap  bifurcation  analysis  the  most  interesting  findings 
are  the  results  associated  with  the  nonlinear  behavior  around  400,000  ft 
starting  altitude.  Of  note  in  this  analysis  is  the  very  marginal 
instability  or  stability  that  may  exist  in  this  region.  In  all  cases  the 
real  part  of  the  eigenvalues  are  very  near  the  imaginary  axis  when  in  this 
altitude  region.  In  the  case  of  the  constant  thrust  rocket  starting  at 
400,000  ft,  the  periodic  branch  was  found  to  have  a  region  of  stability. 
This  implies  that  trajectories  within  a  relatively  small  domain  of 
attraction  would  be  drawn  to  the  limit  cycle  therefore  the  vehicle  could 
expect  to  experience  stable  periodic  motion,  on  the  order  of  the  orbital 
period,  with  fairly  significant  amplitude  for  velocity  and  altitude. 
Looking  at  the  behavior  of  the  translational  states  and  given  that  the 
period  of  the  limit  cycle  is  nearly  equal  to  that  of  a  circular  orbit,  it 
seems  likely  that  the  limit  cycles  associated  with  the  "resonance 
altitude"  describe  the  velocity  and  altitude  of  an  elliptical  orbit. 

The  bifurcation  sweeps  using  the  throttle  parameter  showed  for  the 
most  part  that  the  throttle  is  a  very  benign  way  to  control  the  energy  of 
the  vehicle.  Little  was  found  that  was  of  physical  interest  from  the 
stand  point  of  examining  nonlinear  behavior.  The  most  interesting 
behavior  was  found  for  the  air-breathing  case  with  starting  altitudes 
below  400,000  ft.  At  these  "lower"  altitudes  the  higher  atmospheric 
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density  made  the  aerodynamic  forces  more  effective  and  the  equilibrium 
energy  management  is  done  by  changing  velocity  and  altitude  without 
rotating  the  vehicle.  Only  one  physically  realistic  Hopf  point  was  found 
for  these  low  altitude  sweeps.  The  periodic  motion  looked  similar  to  that 
obtained  from  the  body  flap  sweeps  in  terms  of  the  altitude  and  velocity, 
however  the  rotational  states  experience  a  relatively  high  frequency 
sub-oscillations  relative  to  what  was  seen  before.  In  further  comparing 
the  results  to  the  bifurcation  sweeps  using  the  throttle  parameter  and  the 
body  flap  parameter  for  the  constant  thrust  rocket  case  the  pitch  angle 
and  the  variation  of  the  lift  coefficient  play  key  roles  in  the  dynamics 
of  the  system.  It  seems  without  the  pitch  angle  providing  the  impetus  for 
instability  the  Hopf  bifurcation  at  high  altitudes  where  aerodynamic  pitch 
damping  is  lost  is  not  seen,  which  is  not  altogether  surprising. 

Augmenting  the  vehicle  model  to  obtain  system  stability  for  cases 
where  the  height  mode  is  unstable  seems  to  be  intractable  without  changing 
the  model  to  allow  for  commanded  attitude  input  and  the  ability  to 
generate  or  obtain  the  measurements  of  states  or  some  value  associated 
with  a  state  that  will  allow  for  the  minimization  of  some  error. 

Certainly  the  rich  dynamics  associated  with  nonlinear  phenomena  has 
been  demonstrated  by  the  resulting  complex  behavior  present  even  in  this 
simple  example.  This  work  stands  in  contrast  to  those  in  previous  studies 
who  claimed  to  have  explored  the  nonlinear  nature  of  the  longitudinal 
dynamics  of  a  powered  lifting  hypersonic  vehicle  by  simply  including 
second  order  terms  in  the  Taylor  series  expansions  for  small  perturbation 
analysis.  While  no  great  departures  of  physical  significance  were  found 
in  this  study  from  that  which  was  previously  obtained  by  others,  this  work 
does  display  many  of  the  major  findings  of  previous  works  and  adds  insight 
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to  the  expected  local  behavior  of  trajectories. 

Perhaps  most  useful  of  all,  is  this  work  displays  the  tremendous 
utility  and  encompassing  nature  of  bifurcation  analysis  as  well  as  the 
ease  with  which  it  can  be  applied  to  problems  of  this  sort. 

Recommendations 

This  work  really  stands  as  a  first  step.  It  opens  the  door  to  a 
variety  of  areas  for  further  investigation.  Several  of  these  are  briefly 
mentioned  below. 

1.  Add  the  lateral  equations  of  motion  and  study  the  dynamics  of  a 
powered  lifting  hypersonic  vehicle  flying  a  minor  circle.  This  would 
provide  significant  coupling  of  the  longitudinal  and  lateral  dynamics  and 
should  make  for  some  interesting  behavior. 

2.  Define  the  aerodynamic  forces  and  moment  coefficients  in  nonlinear 
terms  as  found  in  references  such  as  Etkin’s  text  (5:199,393), 

3.  Include  the  rotation  of  the  Earth  in  the  equations  of  motion. 

4.  Increase  the  accuracy  of  the  thrust  laws  to  reflect  more  up-to-date 
propulsion  concepts  such  as  ramjets  and  scramjets.  This  would  bring  Mach 
number  and  additional  altitude  dependencies. 

5.  Develop  higher  order  control  systems  to  stabilize  the  height  mode 
without  destabilizing  the  phugoid  mode.  This  will  require  the  addition  of 
states  to  the  model  to  allow  for  commanded  attitudes  and  feedback  to 
controls  with  direct  influence  over  altitude  and  velocity.  Recall  for 
this  study  the  simple  feedback  of  available  states  to  the  body  flap  proved 
worthless  for  stabilizing  the  height  mode  and  of  little  value  in 
maintaining  phugoid  stability  much  beyond  an  original  Hopf  point. 
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Appendix  A:  FORTRAN  Listing  of  the  User  Supplied  Subroutines  for  AUTO 


CURRENT  AS  OF  23  Nov  1990 


SUBROUTINE  FUNC ( NDIM , U, ICP , PAR , I JAC , F , DFDU , DFDP ) 


This  subroutine  evaluates  stationary  solutions  from  the  equations  of 
motion  for  a  powered  lifting  aerospace  vehicle  flying  along  a  great 
circle  about  a  nonrotating  spherical  Earth. 

Input  parameters  : 


NDIM 

U 


PAR 
PAR ( 1 ) 
ICP 
I  JAC 


Dimension  of  U  and  F. 

State  Vector  containing  U. 

V/VO  the  velocity  along  the  flight  path  divided 
by  a  constant  ( nondimensional ) 

Angle  of  Attack  [radians] 

Pitch  Rate  of  the  Body  relative 
to  the  Earth  [rad/sec] 

Pitch  Angle  [radians] 

Radius  from  Earth's  Center  divided  by  a 
constant  (nondimensional) 

Array  of  parameters  in  the  differential  equations, 
df  Body  Flap  deflection  [degrees] 

PAR( ICP ( 1 ) )  is  the  initial  'free'  parameter. 

*1  if  the  Jacobians  DFDU  and  DFDP  are  to  be  returned 
-0  if  only  F(U,PAR)  is  to  be  returned  in  this  call. 


U(1)« 

V/VO 

by  a 

U(  2  )«* 

alpha 

U(3)~ 

q 

U(  4  )- 

theta 

U(5)« 

r/RO 

i 


Values  to  be  returned  : 

F  -  F ( u , par )  the  right  hand  side  of  the  ODE. 

DFDU  -  The  derivative  (Jacobian)  with  respect  to  U. 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  L, K0 , Ky, mslug , Lift 

DIMENSION  U ( NDIM ) , PAR (20) 

DIMENSION  F(NDIM) , DFDU( NDIM , NDIM ) , DFDP ( NDIM, 20 ) 

COMMON  /CNST/  alt , L, Ky , K0 f re , R0 , VO , gs , IRSTST, ITEST 
COMMON  /FUNVL/  rho, g , Cd , Cl , CM, TPM , DR0DU5 ,DGDU5 , S , mslug 
COMMON  /AERO/  CdO , Cda , CIO , Cla , CMO , CMa , CMdf , CMq 
COMMON  /CFORC/  orbper , Drag , Li f t , altw,W, Fc, Tx , Ty 

.  Set  flag  for  Subroutine  Const  to  use  current  values 

.  of  the  States  rather  than  initial  values. 

ITEST-0 

.  Call  Subroutine  CONST  to  obtain  the  States,  plus  the 

.  necessary  constants  and  functional  values 

CALL  CONST(U) 

************************************************************************** 

**************  System  of  5  Nonlinear  Equations  of  Motion  ***************** 
************************************************************************** 

.  dv/dt  SCALED  ie  U(l)-  V/VO  (Note  U(l)  is  nondimensional) 

.  NOTE:  TPM  -  T/m 
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F(  1 )  ■={  ( TPM )  *DCOS  ( U { 2 ) )  -  ( 0 . 5D0*rho*S*Cd*  ( V0*U(  1 ) )  **2  )/mslug 
&  -g*(DSIN<U(4) }*DC0S(0(2) )-DC0S(U{4) )*DSIN(U(2) ) )  )  /  VO 

..  d( alpha )/dt 

F{ 2 )=U( 3 ) + (  g/(U(l)*V0)  ) * ( DCOS ( U( 4 ) ) *DCOS ( U( 2  ; ) +DSIN( U( 4 ) ) 

&  *DSIN(U( 2 ) ) )-( TPM/(U( 1 ) *V0 ) ) *DSIN(U( 2 ) ' - 

&  ( 0 . 5D0*rho*S*Cl*U( 1 ) *V0/mslug ) 

. .  dq/dt 

F(  3  )«=(  (  0 . 5D0*rho*S*L*CM*  (U(  1 )  *  VO )  *  *  2 )  /  (mslug*Ky**2 )  ) 

&  -(  ( 1 . 5D0 ) * (  g/(U(5)*R0)  ) * ( KO ) *DSIN( 2 . 0D0*U( 4 ) )  ) 

..  d( theta )/dt 

F(4)«U(3)+(  (  (U(1)*V0)/(U(5)*R0)  )*<  DCOS(U( 4 ) )*DC0S(U(2) )+ 
&DSIN(U(4) )*DSIN(U(2) )  )  ) 

. .  dr/dt 

F(5)=( (U{1)*V0)/R0)*(DSIN(U(4) )*DC0S(U(2) )- 
5  DCOS (U( 4) )*DSIN(U(2) ) ) 

IF( I JAC.EQ.O) RETURN 


RETURN 

END 

SUBROUTINE  STPNT ( NDIK , U , PAR ) 


..  In  this  subroutine  the  steady  state  starting  point  must  be  defined. 

..  (Used  when  not  restarting  from  a  previously  computed  solution). 

..  The  problem  parameters  (PAR)  may  be  initialized  here  or  else  in  INIT. 

NDIM  -  Dimension  of  the  system  of  equations. 

U  -  Vector  of  dimension  NDIM. 

Upon  return  U  should  contain  a  steady  state  solution 
corresponding  to  the  values  assigned  to  PAR. 

PAR  -  Array  of  parameters  in  the  differential  equations. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  L, KO , Ky ,mslug 

DIMENSION  U(NDIM) ,PAR(20) 

COMMON  /CNST/  alt , L, Ky , KO , re , RO , VO , gs , IRSTST, ITEST 
COMMON  /FUNVL/  rho , g , Cd , Cl , CM, TPM , DR0DU5 , DGDU5 , S , mslug 
COMMON  /AERO/  CdO , Cda , CIO , Cla , CMO , CMa , CMdf , CMq 

..  Initialize  the  problem  parameters. 

PAR(1)-0,0D0 

write(*,*)  'Enter  initial  par(2)«Kq  and  par(3)~W/S 
&  par(4)»  Trho,  and  par (S)-throttle' 
read ( * , 5 )PAR( 2 ) 
read( * , 5 ) PAR( 3 ) 
read( * , 5 )PAR( 4 ) 
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read( * , 5 )PAR( 5 ) 
5  format(dl5.6) 

U(  2  )«=0 . 0D0 
•U( 4 )=0 . 0D0 
ITEST=10 
CALL  CONST(U) 

C 

RETURN 

END 

C 

SUBROUTINE  INIT 


.  In  this  subroutine  the  user  should  set  those  constants  that  require 

.  values  that  differ  from  the  default  values  assigned  in  DFINIT. 

.  (See  the  main  documentation  for  the  default  assignments). 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

COMMON  /BLBCN/  NDIM, IPS , IRS , ILP, ICP ( 20 ) , PAR( 20 ) 

COMMON  /BLCDE/  NTST,NCOL, I AD, ISP, ISW, IPLT,NBC,NINT 
COMMON  /BLDLS/  DS , DSMIN , DSMAX, IADS 
COMMON  /BLLIM/  NMX, NUZR, RLO , RLl , AO, Al 
COMMON  /BLMAX/  NPR, MXBF , I ID, ITMX, ITNW, NWTN , JAC 

★  ***************************a**;*A******************ifrdrri:*ifc****ifc******* 

C* ******************  KEAD  AUTO  PARAMETERS  *************************** 
C******************************************************************** 

C 

OPEN(UNIT=27, FILE- 'DS. DAT' , STATUS* 'OLD' ) 

REWIND  (27) 

C .  ITEMS  IN  COMMON  BLBCN  -  BASIC  CONSTANTS 

READ( 27 , * )  NDIM 
READ(27,*)  IPS 
READ( 27 ,  *  )  IRS 
READ( 27, * )  ILP 
READ( 27, * )  ICP(l) 

write(*,*)  'Which  parameter  to  vary?(PAR(2 )«Kq  PAR(3)»W/S 
£.  PAR(4)-Trho,  PAR( 5 )-Throttle  setting)' 


read(*,5)  I 
5  format(Il) 

ICP(2)=I 

C .  ITEMS  IN  COMMON  BLCDE  -  DISCRETIZATION  CONSTANTS 


READ( 27,  * )  NTST 
READ( 27 ,  * )  NCOL 
READ( 27 , * )  IAD 
READ( 27 , * )  ISP 
READ( 27 , * )  ISW 
READ(27,*)  IPLT 

C .  ITEMS  IN  COMMON  BLDLS  -  STEPSIZE  ALONG  SOLN  BRANCHES 

READ( 27 , *  )  DS 
READ( 27 , * )  DSMIN 
READ(27,*)  DSMAX 
READ (27,*)  IADS 

C .  ITEMS  IN  COMMON  BLLIM  -  LIMITS 

READ (27,*)  NMX 
READ( 27 , * )  NUZR 
READ( 27 , * )  RLO 
READ( 27 , * )  RLl 
READ( 27 , * )  AO 
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READ( 27 , * )  Al 

C .  ITEMS  IN  COMMON  BLMAX  -  MAXIMA 

READ (27,*)  NPR 
READ (27,*)  MXBF 
READ( 27 , * )  IID 
READ (27,*)  ITMX 
READ( 27 , * )  ITNW 
READ (27,*)  NWTN 
READ (27,*)  JAC 
CLOSE  (27) 

C 

RETURN 

END 

C 

FUNCTION  USZR( I ,NUZR,PAR) 


.  This  subroutine  can  be  used  to  obtain  plotting  and  restart  data 

.  at  certain  values  of  free  parameters. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  L, Ky, KO ,msiug 

DIMENSION  U ( 5 ) , PAR (20) 

COMMON  /CNST/  alt , L , Ky , KO , re , R0 , VO , gs , IRSTST, ITEST 
COMMON  /FNVL/  rho,g,Cd,Cl,CM,TPM,DRODU5,DGDU5,S,mslug 

.  Initially,  for  the  steady  state  analysis,  set  NUZR-0  in  INIT. 

.  Then  the  functions  specified  below  will  be  ignored. 

.  When  computing  the  branch  of  periodic  solutions,  set  NUZR-4  in  INIT. 

.  Output  will  then  be  written  in  unit  8  for  the  values 

.  of  PAR ( * )  specified  below. 

.  Note  that  PAR(ll)  is  normally  reserved.  It  is  used  by  AUTO  to  keep 

. .  track  of  the  period  (See  main  documentation). 

GOTO( 1 , 2 , 3 , 4 ) I 

1  USZR*=PAR(  11 )  -  10.0 

RETURN 

2  USZR«PAR( 11 )  -  14.0 

RETURN 

3  USZR=PAR( 11 )  -  20.0 

RETURN 

4  USZR-PAR(ll)  -  30.0 

RETURN 

END 

SUBROUTINE  CONST(U) 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  L,K0,Ky,Lift,mdot,m,mslug 
C 

DIMENSION  U(NDIM) 

COMMON  /BLBCN/  NDIM, IPS , IRS , ILP, ICP( 20 ) , PAR ( 20 ) 
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COMMON  /CNST/  alt,L,Ky,K0, re,R0,V0,gs, IRSTST,ITEST 
COMMON  /FUNVL/  rho, g, Cd , Cl , CM, TPM, DRODU5 ,DGDU5 , S ,mslug 
COMMON  /AERO/  CdO , Cda , CIO , Cla , CMO , CMa , CMdf , CMq 
COMMON  /CFORC/  orbper ,Drag , Lif t , altw, W, Fc, Tx, Ty 
C 

C**************  Aerodynamic  and  Geometric  Constants  *********************** 
£******★*******************************************£******■***************** 
Cd0=0 . 0133D0 
Cda«=0 . 4D0 
C10=0 . 05D0 
Cla=0.5D0 
CM0=0.0D0 
CMa=-0 . 0548D0 
CMdf=CMa*l . 5D0 
CMq=-0 . 028D0 
C 


C .  Characteristic  Length  of  Vehicle  -  overall  length  L=50  ft 

C 

L=50 . 0D0 

C . .  Weight,  Mass  and  Area  of  Vehicle 


C 

Wsl=700 . 0D3 


S“Ws 1/PAR ( 3 ) 
mslug=Wsl/32 . 174D0 
C 

C .  Radius  of  gyration  in  pitch  [ft]  -  Ky~2  «  Iy/m 

C 

Ky=25 , 0D0 
C 

C .  K0°( Ix-Iz )/Iy 

C 

K0=-0 . 94D0 
C 

C .  Radius  of  the  Earth  [ft]  (standard  geoid) 

C 

re-2.0903264468D7 

C 

C......  Gravity  at  Earth's  surface  [ft/sec‘2] 

C 

gs=32 . 174D0 
C 


V0«1000.0D0 
R0=1000 . 0D0 

raddeg»2 . 0D0*3 . 14159265359D0/360 . 0D0 
C 

C  ************************************************************************* 

C  ************  read  altitude  and  set  initial  U ( 5 )  ************************* 
c  ************************************************************************* 
c 

X F(XTEST  .GT.  5)  THEN 

OPEN ( UNIT“25 , FILE® 'ALT. 25 ' , STATUS= ' OLD' ) 

REWIND  (25) 

READ( 25, * )  alt 
C 

C .  SET  INITIAL  VALUE  FOR  THE  RADIUS  U(5).  NOTE:  SCALED  DOWN  BY  RO 

C 

U( 5 )«( re+alt )/R0 
C 


*nn  non  n  non  n  non  nnnn  n  n  nnnn  non  nnnnn 


■ 


************************************************************************* 
**************  Density  and  Gravity  variation  with  altitude  *************** 

rt************************************************************************ 


I F (  (U(5)*R0  -  re)  .LT.  6.0D5  )  THEN 
..  Convert  altitude  to  SI  units  {km} 

altsi=(  <  U(5)*R0  -  re  )  /  3.2808399D0  )  /  1000. 0D0 
..  Calculate  the  density  in  SI  units  {Kg/nT3] 

..  Constants  used  in  polynomial  (U.S.  Standard  Atmosphere  Supp.  1966) 

co=o.ioooooooood+oi 

C1=0.3393495800D-1 
C2=-0 . 3433553057D-2 
C3=0.5497466428D-3 
C4  — 0.3228358326D-4 
C5=0.1106617734D-5 
C6=-0 . 2291755793D-7 
C7=0.2902146443D-9 
C8—0 . 2230070938D-11 
C9=0 . 1010575266D-13 
CIO— 0.2482089627D-16 
Cl 1=0 . 2548769715D-19 

eqnl=CO+(Cl*altsi )+(C2*altsi**2 )+(C3*altsi**3 )+{ C4*altsi**4 )+ 
&(C5*altsi**5)+(C6*altsi**6)+(C7*altsi**7)+(C8*altsi**8)+ 

&( C9*altsi**9 )  +  (C10*altsi**10 )+(Cll*altsi**ll ) 

rhoSI=l .2250D0/eqnl**4 

..  Convert  from  Kg/nT3  to  Lbra/Ft“3 
..  (0.062427961  ( lbm/f tA3 )/( Kg/jn*3 )  ) 

rho=rhoSI*0 . 062427961D0 

..  Convert  to  [Slugs/ft'3] 

rho=rho/32 .174D0 

ELSE 

..  Equation  for  rho  if  altitude  is  greater  than  600,000  ft 
rho=  2. 16871253724 D- 10  * 

&  DEXP{-8 . 89837671693 D-6  *  (U(5)*R0  -  re)  ) 


END  IF 

.  Calculate  the  gravitational  acceleration 

g-  gs  *  (  re  /  (U(5)*R0)  )**2 

.  Calculate  the  Atmospheric  Pressure 
p0=2116.22D0 

*  pi  —  5.850746831820396D-05 

*  p2=9. 792179784 448163D- 01 

*  p3=9. 875326461241002 D- 05 


81 


on  ono  *-  *  *-  o  o  o  n  non  o  o  00000*000*00  00000 


p4=-6.044333173347913D-06 
p5=3.408653276857509D-09 
p6=-8 .934489792146698D-07 
alt=(U{ 5 ) *R0  -  re) 

Pa=pO*DEXP(pl*alt**p2)+p3*DEXP(p4*alt)+p5*DEXP(p6*alt) 

****************************************£****★**&*********************★*** 

*******  Thrust  Equations  and  Aerodynamic  Coefficients  ******************** 
************************************************************************** 

Cd=CdO+Cda*U( 2 ) **2 
Cl=ClO+Cla*U( 2 ) 

.  Exhaust  Nozzle  Area  [ft"2] 

An=40 . 0D0 

.  Exhaust  Velocity  (Vexh)  *  500  ft/sec 

Vexh=500.0D0 

.  Rocket  with  constant  Thrust  at  reference  altitude  and  velocity. 

. Note:  W/S-30  [lb/ft*2]  at  sea  level,  therefore 

.  U(l)initial  -  sqrt{  (g  r)  /  (1  +  ( rho  r  Cl  S)/2m))  is 

IF  (ITEST  .gt.  5)  THEN 

U(l)*=  DSQRT(  ( g*U( 5 ) *R0 )/(  1 . 0D0+( rho*S*U( 5 ) *R0*C1 )/ 

&  (2.0D0*mslug)  )  ) 

U(3)  -  -U(1)/(U(5)*R0) 

.  Thrust  Constant 

T0=(  0.5D0*rho*S*Cd*U(l)**2  )/<  rho**PAR(’4 ) ) 
write(*,*)  ' TO- ' , TO 

.  Mass  flow  rate  in  (slugs/sec] 

.  Assume  the  Pressure  is  expanded  to  the  starting  altitude  value 

Pe*»Pa 

mdot=(  (0.5D0*rho*(U(l)**2)*Cd*S)  -  (Pe-Pa)*An  )  /  Vexh 
write(*,*)  'MOOT*' , indot 

.  SCALE  DOWN  U(l)  and  nondimensionalize 

U(1)-U(1)/V0 

.  PRINT  RESTART  DATA  FILE 

OPEN (27,  FILE** 'REF. DAT'  .STATUS** 'NEW' ) 

WRITE ( 27,10)  alt, TO 

1C  FORMAT(40x, 'REFERENCE  VALUES' ,/,6X, 'ALT  [FT]', 

&  14X, 'TO  [ f t‘4/sec*2 ] ' ,/, 2 ( IX, E15. 8 , 4X) ) 

CLOSE( 27 ) 

OPEN(  40, FILE- 'THRUST. DAT'  , STATUS** ' NEW'  ) 

OPEN( 41 , FILE- ' COEFF.DAT' , STATUS- ' NEW' ) 

WRITE( 40,21) 

WRITE( 41,22) 

IRSTST-10 
try=0. OdO 

ELSEIF  ( IRSTST  .LT.  5)  THEN 
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OPEN(27,FILE='REF.DAT' ,STATUS-'OLD* ) 

REWIND{ 27 ) 

READ( 27 , 20 )  x jnkl ,T0 

OPEN ( 40, FILE®' THRUST. DAT' , STATUS® 'NEW' ) 

OPEN( 41 , FILE= ' COEFF.DAT' , STATUS* 'NEW' ) 

WRITE(40,21) 

WRITE( 41,22) 

20  FORMAT (/,/, 2(1X,E15.8,4X) ) 

21  FORMAT ( 6X, 'H' ,15X, 'Tx' ,15x, 'D' ,15x, 'Ty' ,15x, 'L' ,15x, 'W' , 

&  15x, 'Fc' ,14x, 'Orb  Per') 

22  FORMAT(8X, 'H' ,15x, 'Cl' ,15x, 'Cd' ,15x,'CM' ,15x, 'V' ,15x, 

&  'alpha' ,15x, 'Rad' ) 

IRSTST=10 
END  IF 
C 

Thrust  =  (T0*rho**PAR( 4 ) )  *  PAR(5) 
orbper®  ( 6 . 2831853072D0  *  U(  5 ) *R0 )/DSQRT( g*U( 5 ) *R0 ) 

Drag®  0.5D0*rho*S*Cd*(U(l)*V0)**2 
Lif t=  0.5D0*rho*S*Cl*(U(l)*V0)**2 
altw=  U(5)*R0  -  re 
w  =  mslug*g 

Fc  =  mslug* (  (  (U(1)*V0)**2  )  /  (U(5)*R0)  ) 

Tx=Thrust*DCOS(U(2 ) ) 

Ty=Thrust*DSIN(U(2) ) 

TPM=Thrust/raslug 

if  (try  .It.  l.OdO)  write(*,*)  'TPM=',TPM 
try«=5.0d0 

.  PAR ( 1 )  -  body  flap  deflection 

.  qO  -  Pitch  rate  due  to  spherical  Earth  in  body  axis  system 

.  Note  all  angles  are  in  radians  EXCEPT  the  body  flap  deflection 

.  (  df=PAR( 1 )  )  which  is  in  degrees,  therefore  PAR(l)  is  multiplied 

.  by  2*pi/360  [ rad/deg) 

qO  =  - ( U( 1 )  *  VO )/(U( 5 )  *  R0) 

U0  -  DSQRT(  g*U(5)*R0/((rho*U(5)*R0*Cl/2.0D0*mslug)+1.0D0)) 
df  -  PAR ( 1 )  -  PAR(2)*(U(3)-qO) 

CM  «*  CM0+CMa*U(2)+CMq*(U(3)-q0)+CMdf*raddeg*df 

RETURN 
END 

subroutine  BCND 
return 
end 

subroutine  ICND 
return 
end 

subroutine  FOPT 
return 
end 
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Appendix  B:  Standard  Atmospheric  Approximations  for  Density  and  Pressure 


The  value  of  density  for  the  Standard  Atmosphere  (14)  is  calculated  using 
a  different  equation  over  two  altitude  regions.  The  lirst  altitude  range 
is  from  0  to  600,000  ft.  The  density-altitude  approximation  for  this 
range  was  obtained  directly  from  the  work  by  Vihn  and  Dobrzelecki  (15:25) 
and  provides  values  of  density  accurate  to  within  5%  of  the  Standard 
Atmosphere  for  altitudes  ranging  from  0  to  200  Km  (0  to  656,000  ft) 
(15:25).  This  equation  is  an  inverse  polynomial  relationship  given  by: 


_ P£J _ 

[A0  +  \Z  +  ...  +  A11Z11]  * 


(41) 


where: 


P 


Psl 


Z 


density  [kg  m’3] 

sea  level  density  [kg  m’3] 

altitude  above  standard  geoid  (6371.315  km) 

(note  this  is  the  average  radius  of  the  Earth  at  the 
equator,  which  is  different  for  reference  15) 
Coefficients  (j=l-ll)  [km'3] 


Table  2 

Coefficients  for  Density  Polynomial 


0 

0.1000000000 

E  01 

1 

0.3383495800 

E-01 

2 

-0.3433553057 

E-02 

3 

0.5497466428 

E-03 

4 

-0.3228258326 

E-04 

5 

0.1106607734 

E-05 

6 

-0.2291755793 

E-07 

7 

0.2902146443 

E-09 

8 

-0.2230070938 

E-10 

9 

0.1010575266 

E-13 

10 

-0.2482089627 

E-16 

11 

0.2548769715 

E-19 

(15:26) 


r 


;  Note  the  density  is  computed  in  [kg  m"3]  and  is  then  converted  to  [slugs 

ft'3]  using  1.9403232735  E-03  [{kg  m'3)/ (slugs  ft"3)]. 

For  altitudes  above  600,000  ft  the  following  exponential  relation  is  used: 

p  =  2. 16871253724E-10  exp  (-8. 898376716935-06  Z)  (42) 


where:  p  =  density  [slugs  ft"3] 

Z  =  altitude  [ft] 

Note  this  equation  yields  the  value  of  p  in  [slugs  ft"3]  , 

Figure  44  shows  the  calculated  density  relative  to  the  Standard 
Atmospheric  data  from  reference  14. 

For  the  pressure  altitude  relation  an  exponential  relationship  was 
developed  using  the  nonlinear  fitting  routine  on  MATLAB  from  MathWorks. 


P  =  PO  ExpiPl  ZP2)  +  P3  Exp(P4  Z)  +  P5  Exp  { P6  Z)  (43) 


where:  P 

P0 
PI 
P2 
P3 
P4 
P5 
P6 
Z 


atmospheric  pressure  [lb  ft  ] 

2116.22  [lb  ft'2] 

-5.850746831820396  E-05 
9.792179784448163  E-01 
9.875326461241002  E-05  [lb  ft"2] 
-6.044333173347913  E-06 
3.408653276857509  E-09  {lb  ft'2] 
-8.934489792146698  E-07 

equatorial  altitude  above  the  Standard  geode  [ft] 


Figure  45  shows  the  quality  of  the  pressure  fit  to  the  Standard 
Atmospheric  data  in  reference  14. 
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Pressure  [psf ]  c  0ensity  {3(ug/cu  ft] 
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